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Abstract 




The three-dimensional, time-dependent (incompressible) vorticity 
equations have been used to simulate numerically the decay of isotropic 
box turbulence and time-developing mixing layers. The vorticity equa- 
tions are spatially filtered to define the large-scale turbulence field, 
and the subgrid scale turbulence is modeled. A general method has been 
developed to show numerical conservation of momentum, vorticity, and 
energy that is much simpler than previous methods and is widely appli- 
cable. The terms that arise from filtering the equations have been 
treated (for both periodic boundary conditions and no-stress boundary 
conditions) in a fast and accurate way by using fast Fourier transforms. 
Use of vorticity as the principal variable is shown to produce results 
equivalent to those obtained by use of the primitive variable equations. 

A new subgrid scale model is used in conjunction with the vorticity 
equations and is shown to produce results that compare well with the ex- 
perimental results. The new model offers advantages both in computational 
speed and in storage. 

The vortex-pairing mechanism, observed in the spatially developing 
counterpart of the time-developing mixing layer, has been simulated nu- 
merically. It is interesting to note that with simply two vortices pair- 
ing, self-similar mean velocity and mean turbulence intensity profiles 
are obtained. The vortex-pairing mechanism is shown to be persistent 
even with the presence of large-amplitude, three-dimensional background 
turbulence. A number of different initial fields have been studied. The 
presence of large organized structures, in the initial conditions, is 
shown to be essential in order to predict growth rates of the mixing lay- 
ers comparable to those observed experimentally. The rate of growth is 
found to be very dependent on the initial field, a fact also observed 
experimentally. 
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Chapter 1 
INTRODUCTION 

1.1 Background 

fi . ’ .'•'•• '■■'.••• ' • : '• ••.' 

Turbulent flows have been the subject of experimental and theoretical 
investigations since the last century. Despite the formidable amount of 
effort invested in this field, our ability to predict flows of technical 
importance remains severely limited. The major difficulty encountered by 
the theoretical investigations arises from the nonlinear character of the 
equations of motion. Statistical averages of the equations of motion give 
rise to the so-called Reynolds stresses. The equations for the Reynolds 
stresses in turn give rise to higher-order statistical quantities, and so 
on. The usual approach to computing turbulent flows is to model- the terms 
that arise from the nonlinear character of the equations of motion. This 
approach usually requires a great deal of experimental information. 

We know the underlying physical principles of most fluid flows, and 
the quantities of interest are completely determined by known equations. 
With the introduction of large computers, three-dimensional, time-dependent 
computation of turbulent flows has become possible. However, in order to 
resolve all the scales of motion even in the simplest turbulent flow, 
namely, the isotropic homogeneous case, Kwak et al. (1975) estimated the 
number of mesh points needed in any given direction to be 



( 1 . 1 ) 


where'. 

Rt “ (qJf/v) , 

V = kinematic viscosity, ; 

X - length scale of large eddies, and 
q = r.m.s. velocity. 

Equation (1.1) shows that one can do a full simulation only at very 
low Reynolds number. Indeed, Clark et al. (1977), using a 64 x 64 x 64 
mesh system, were able to solve the isotropic homogeneous turbulence 


problem for - qX/v = 38.1, where X is the Taylor microscale. Their 
predicted results compared well with the experimental results. However, 
turbulent flows of technical importance have much higher Reynolds num- 
bers, and all the scales of motion cannot be resolved for these flows. 

One of the more promising approaches to solving turbulence problems 
is "large-eddy simulation". In large-eddy simulation, one calculates 
the large-scale turbulent motions with a relatively coarse time-dependent, 
three-dimensional computation that uses some sort of model (the "subgrid 
scale model") for the small scales. The basic motivations for this ap- 
proach are twofold. First, experimental observations of turbulent flows 
show that the large turbulent structures differ markedly from one flow 
type to another (e.g., jet vs. boundary layer), but the small-scale tur- 
bulent structures are quite similar. Thus, while there is little hope 
of concocting a "universal" model for the large structures, it may be 
possible to do so for the small-scale motions. Second, as computer 
capabilities grow, our capability of resolving smaller scales Will grow 
and the effects of the subgrid scale model will diminish. Thus, while we 
are limited to simple flows with the present computer capabilities, large- 
eddy simulation is a tool that may be used on future generation computers. 

Kwak et al. (1975) and Shaanan et al, (1975) have shown that homoge- 
neous turbulent flows can be simulated reasonably well with a relatively 
small number of mesh points (16 x 16 x 16) . Orszag and Pao (1974) , using a 
32 x 32 x 32 mesh system, predicted the momentumless wake of a self- 
propelled body. Deardorff (1970) and Schuman (1973) computed the central 
region of a plane channel flow using the large-eddy simulation approach. 
While Deardorff and Schumann did not handle the wall (no slip) problem, 

Moin et al. (1978) have solved the channel flow problem, including the 
laminar sublayer. In this work we shall study the time-developing, two- 
stream mixing '.layer:.: . 

Previous works on prediction of the two-stream mixing layer have con- 
centrated on the initial stages (roll-up) of the development of the layer. 
Patnaik et al. (1976), starting with an initial distribution that is an 
unstable eigensolution of the Taylor-Goldstein equation, predicted the 
two-dimensional roll-up of a stably stratified horizontal mixing layer. 
Another method that has been used to compute the mixing layer in two 

2 



dimensions is the vortex-tracing method used by Ashurst (1977). This 
method suffers from high computational costs and ad-hoc assumptions con- 
cerning the effects of viscosity. The high computational cost of the 
vortex-tracing method can be reduced by using the vortex- in-cell method 
(Wang, 1977). These works have treated two-dimensional cases, but the 
mixing layer exhibits three-dimensionality. This is apparent from the 
shadowgraph pictures of Brown and Roshko (1974) and the spanwise velocity 
fluctuation measurements of Spencer and Jones (1971). 

1.2 Experimental Background 

The two-dimensional turbulent mixing layer plays an essential role 
in many technological problems . For example, the initial regions of pla- 
nar jets can be approximated as two independent, two-dimensional mixing 
layers. Flow over a backward-facing step (with a large step height) is 
another example of the two-dimensional mixing layer. Many other flow 
situations can be identified with the mixing layer. In combustion pro- 
cesses, fluid mechanics plays a major role in mixing the reactants, and 
better understanding of turbulent: mixing is needed. The mixing layer is 
perhaps the simplest situation in which two flows come into contact; ob- 
viously the ability to analyze simple problems is necessary before one can 
analyze more complicated ones. 

In 1947, Liepmann and Laufer studied the mixing layer and established 
the general features of the flow. However, the fundamental understanding 
of the structure of the flow is still far from complete, and many contro- 
versial questions need to be answered. We shall address some of these 
questions. The reader is referred to Mur thy (1975) for an extensive re- 
view and interpretation of the available literature on the mixing layer. 
With the advancement of the techniques of hot-wire anemometry, Wygnanski 
and Fiedler (1970) attempted to reproduce Liepmann and Laufer data and 
extend it to include other measurements. However , differences in inten- 
sity levels and rate of growth of the layer emerged. These differences 
were attributed to the presence of a trip wire in the Wygnansky and Fied- 
ler experiment that was not used by Liepmann and Laufer. Batt (1975) 
studied both configurations and showed that the differences are due to the 
tripping of the layer. Foss (1977) investigated the effects of the 
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laminar/turbulent boundary layer states on the development of a plane 
mixing layer. He found that the development of the layer is dependent on 
the initial conditions (che status of the boundary layers before the two 
streams merge) . Figs. 1.1 and 1.2 show r.m.s. fluctuations of the stream- 
wise velocity and the mean velocity profiles obtained by Foss. These 
figures show that different self-similar stages are obtained for different 
initial conditions. Foss argues that this is due to the sensitivity of 
choosing the virtual origin of the mixing layer (x q ) and that the charac- 
ter of the (inif: •!) disturbance, not its amplitude, is responsible for 
the substantial effect on the virtual origin. More recently, Oster et al. 
(1977) showed that by oscillating the initial conditions of the mixing layer 
they can more than double the growth rate of the layer. The effect depends 
on the frequency and amplitude of the oscillations introduced. These ex- 
perimental results show that the "universality" of the self-similar stage 

.. . 6 ; 

of the mixing layer is in doubt, at least up to Re = 1.5 x 10 . Fiedler 
and Thies (1977) showed that the two-dimensional shear layer only slowly 
reaches a self-similar state and that every disturbance is of long influ- 
ence, Table 1.1 shows tabulated results extracted from the Fiedler and 
Thies paper, and it can be clearly seen that different experiments predict 
different growth rate of the layer. 

Winant and Browand (1974), using dye visualization in a mixing layer, 
observed that initially the fluid rolls up into discrete, two-dimensional 
vortical structures. These structures then interact by rolling around 
each other to form a single larger structure. This pairing process con- 
trols the growth of their mixing layer. Brown and Roshko (1974) also ob- 
served the amalgamation process at Reynolds number 2.5 x10 . Chandrsuda 
and Bradshaw (1975) argue that the two-dimensional, large-eddy structure 
observed by Brown and Roshko is unlikely to survive indefinitely if the 
ambient entrained fluid is weakly turbulent. They advance the argument: 

"It is probable that if the Brown and Roshko type of orderly structure is 
once formed it can last for a large number of characteristic wavelengths — 
that is, up to high Reynolds numbers based on streamwise distance — but 
not indefinitely. The question can be settled only by measurements in a 
two-stream mixing layer at a much higher Reynolds number than was used by 
Brown and Roshko." Dimotakis and Brown (1976) showed the existence of 
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large structures at Reynolds number - 3 * IQ and attributed the growth 
of the mixing layer to both pairing and "tearing* 1 . Tearing is described 
in their paper as an event where "a large structure will occasionally 
find itself in the vicinity of another, or in between two others, in whose 
straining field it disintegrates.*' The tearing process was first advanced 
by Moore and Saffman (1975) on the basis of exact solutions for uniform 
vortices in straining fields. 

1.3 Motivation and Objectives 

In many flows of practical interest there are interactions between 
irrotational regions and turbulent regions. Examples of such flows are 
the shear layer, turbulent jets, and turbulent boundary layers with irro- 
tational free stream flow. In such flows, the regions are separated by 
a very thin super layer across which there is normally a jump in the vor- 
ticity components parallel to the layer. The dynamical equations for the 
vorticity seem to be suited to simulate such flows, since the vorticity is 
identically zero in the irrotational region. However, previous workers 
used the dynamical equations in the primitive variables (velocity, pres- 
sure) and there has been doubt (Orszag and Israeli, 1974) that the vor- 
ticity equations could be used to solve turbulent flow problems. Our 
objectives were therefore as follows: 

• To explore the feasibility of using the vorticity equation to simu- 
late turbulent flows. 

• To find a subgrid scale model appropriate to the vorticity equations 
and to determine any constants in this model. 

• To simulate a turbulent flow with interactions between turbulent re- 
gions and non-turbulent irrotational regions; we chose the mixing 

V layer. 

In order to use the three-dimensional, time- dependent vorticity equa- 
tions, we need to develop a numerical approximation based on these equa- 
tions that conserves mass, momentum, vorticity, and energy . We also need 
to assess numerical finite-difference methods and, in particular, the 
fourth-order and pseudo-spectral methods. 


1 . 4 Overview 


The equations of motion of the large eddies are derived by averaging 
(filtering) the vorticity equations in space. In Chapter 2, we describe 
the approach to solving turbulent flow problems that is called large-eddy 
simulation. We show that the use of a filter that is smooth in the real 
space is required to handle rotational- irrotational regions. We present 
a new subgrid scale model to be used in conjunction with the vorticity 
equations that is much simpler and faster than the one that would be ob- 
tained from the more commonly used Smagorinsky model. 

In Chapter 3, we describe the numerical methods used in this work, 
briefly discussing the fourth-order and pseudo-spectral approximations and 
numerical filtering. We develop a numerical approximation to the vorticity 
equation that conserves mass, momentum, vorticity, and energy, and a method 
of deriving conservation properties that is much simpler than previous 
methods and is widely applicable. We present a new treatment of the fil- 
tered convective and stretching terms that is more accurate and faster 
than previously used methods. 

In Chapter 4, the isotropic homogeneous turbulence problem is solved 
using both fourth-order differencing and the pseudo-spectral approximation. 
The numerical approximations to the partial derivatives of the subgrid scale 
model are discussed. We show that the use of the vorticity equation to 
solve turbulent flow problems is feasible and that the new model produces 
results equivalent to those produced by previously established models. 

In Chapter 5, we discuss the two-dimensional computation of a mixing 
layer. An array of vortices is perturbed, and the momentum thickness 
growth rate is discussed as a function of the perturbation. It is inter- 
esting to note that self-similar, mean velocity and turbulence intensity 
profiles are obtained with vortex pairing. 

In Chapter 6 a three-dimensional computation of a turbulent mixing 
layer is studied. It is found that the presence of large structures in 
the initial conditions is essential for the successful prediction of tur- 
bulent mixing layers. Our studies of different initial conditions pro- 
duce different growth rates of the layer — a fact supported experimentally 
Self-similar, mean-velocity profiles are obtained with different flow struc 
tures. However, turbulence intensity profiles show a rapid decay when 
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large turbulent structures are not present. We show that our subgrid 
scale model inhibits the production of turbulent fluctuations when we 
start with random turbulent fluctuations added to a mean velocity pro- 
file, i.e., the model is incapable of handling transitional flows, with 
present computational limitations. 


Chapter 2 


THEORETICAL FOUNDATIONS 

2.1 Definitions of the Large and Subgrid Scales 

In the previous chapter it was shown that, due to computer limita- 
tions, one cannot do a full simulation of the dynamical equations of tur- 
bulent fluid motion except at extremely low Reynolds numbers. We pointed 
out that the large-scale turbulent structures differ markedly from one 
flow to another (e.g., jet vs. boundary layer), while the small-scale 
turbulent structures are quite similar, and that large-eddy simulation 
is a promising approach. 

In the large-eddy simulation approach, the first and most fundamental 
step is defining the large-scale field. A general approach that recognizes 
the continuous nature of the flow variables is the "filter function" ap- 
proach of Leonard (1973). If f is some flow variable, we can decompose 
it as follows: 

■ ^ . f - f + f ' (2.1) 

where f is the large-scale (filtered) component and f 1 is the residual 
field. Leonard defined the filtered field by: 

= f G(x-x') f(x') dx’ ■■■■/!: -V; • (2.2) 

where G(x-x') is the filter function., and the integral is extended over 
the whole flow field. One can think of f as a local spatial-averaged 
field. 

It can be shown that if G is piecewise continuously differentiable 
and G(r) goes to zero as r -> °° and is integrable over an infinite do- 
main, then 


However , in general 


fg 1 fg (2.3b) 

Properties (2.3a) and (2.3b) will be used in deriving the dynamical 
equations of the large scales motion. 


2.2 Dynamical Equations in Vorticity Form 

In Chapter 1, we pointed out that in many flows of practical interest 
there are interactions between irrotational regions and rotational turbu- 
lent regions. Examples of such flows are shear layers, turbulent jets, 
and turbulent boundary layers with irrotational free streams. In such 
flows the regions are separated by a very thin super layer across which 
there is normally a jump in the vorticity parallel to the layer. These 
flows are a challenge to the experimentalist; the difficulties arise from 
the fact that it is hard to determine the region of the flow in which the 
measurements are made. One faces a similar problem in trying to simulate 
such flows numerically. The difficulty arises from the fact discussed 
earlier, that it is impossible to capture all of the scales of motion in 
the turbulent region. The best we can do is to filter the dynamical equa- 
tions to obtain equations that describe the behavior of the large eddies, 
and to model the small scales. Since in the irrotational region the vor- 
ticity is identically zero, the dynamical equations for the vorticity seem 
to be suited to simulate such flows. 

Now let us derive the dynamical equations for the large-scale vortic- 
ity field. For an incompressible fluid with constant viscosity, the equa- 
tions of motion for the primitive variables may be written: 


9u ■ .. 

i 

-T7 £ . U.to. 

at ijk J k 


9x i W 

3u . v 
x 


§ Vi) 


- ve 


0), 


ijk 3x. k 
3 


3x. 


= 0 


(2.4) 

(2.5) 


The vorticity equation is obtained by taking the curl of Eqn. (2.4). 
Operating on it with e . 3/9x gives : 

- pqX ' q ■ /'V 
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V 


( 2 . 6 ) 


9 w, 


£“i + -dr Vi'W ■ '’a— 

3 3 3 

Multiplying Eqn. (2.6) by a filtering function G(x-x/) and inte- 
grating over the whole flow field, we obtain: 


g g ^ ^4 

0), + “ — (u.W. - U.OO.) = V 


r> vju , i ,v \ u , w , u , 

9t i 3x, ] i i 3 

3 


9x. 9x. 
3 3 


(2.7) 


The fact that a finite-difference approximation of Eqn. (2.7) would 
involve approximating higher derivatives of the velocity than would be the 
case with the primitive equations (Orszag and Israeli, 1974) need not worry 
us in this case. Since the equations are filtered, we shall be dealing 
with smooth functions. 

As can be expected, when averaging nonlinear equations, we run into 


the closure problem; i.e. , we need to express the quantities u^Oh and 
u^w_, in terms of u and 0 ). Expanding u and U) as in Eqn. (2.1), 
one obtains 


U.GJ. - U .Cl) , 

3 1 1 3 


u .0) . - u.co. + W. . 
3 1 13 13 


( 2 . 8 ) 


where 


W.. = u.(i)’ + u!tu. - u.u)! - u!u). + u!o)! - u!u>! 

13 3 i ji 13 13 31 13 


(2.9) 


We note that W.. contains subgrid scale quantities and hence must be 
modeled. 


2.3 Subgrid Scale Models 

We first note that the model of W^. should satisfy the following 
necessary conditions': ; 

1. Antisymmetry, since W^. is an antisymmetric tensor and therefore 


W, 


9x , 9x . ij 
i 3 


= 0 


( 2 . 10 ) 


It is important to preserve the antisymmetry property of W^j in 
order to assure 9u>^/9x^ = 0, since the dynamical equations for the 
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vorticity do not contain a pressure-like term which could be used 
to adjust the divergence of the vorticity. 

2, It should vanish in an irrotational region, since W,. vanishes in 

*^*3 

such regions. 

3. It should be an energy sink, since it represents subgrid-scale effect 


2,3.1 Model co - 1 

Previous workers (Kwak et al. , 1975; Shaanan et al . , 1975), working 
With the filtered dynamical equations in the primitive variables* used an 
eddy-viscosity model for their subgrid-scale model. They modeled the term 


t.. - u!u! + u!u. + u!u. - i (u.’u,' + 2u.' u, ) <$.. 

i -3 i 3 3 x i 3 3 k k k W 13 


by setting 


T . . = - 2V„S , , 

13 T 13 


( 2 . 11 ) 


where 


’ij 


( al: u : + 33- u i ) 

N 1 3 


( 2 . 12 ) 


is the strain rate tensor of the filtered field and is an eddy vis- 

cosity associated with the sub grid scale motions, 

Smagor insky (1963) suggested a model for V T 


<«.«* CM S 


(2.13) 


where C is a constant and A is the filter width. We note that in a 
s 

non^ turbulent region this model of may have a non-zero value, and 

hence it may give rise to residual stresses. Since our main objective is 
to handle interactions between a turbulent region and a non- turbulent re- 
gion, this model was rejected for the. present work. 

One way to avoid this difficulty is to relate directly to vor- 

ticity. Previous workers (Kwak et al. , 1975; Donaldson, 1972) used 


V T = (C v A) 2 (Wlih)* 5 (2.14) 


where C v is a constant. Clark et al. (1977) have shown that this model 
is as accurate as Smagorinsky 1 s for homogeneous isotropic turbulence. 

The dynamical equations for large-scale vorticity field could have 
been derived by taking the curl of the filtered dynamical equations for 
the primitive variables. Hence the curl of Eqn. (2.11) could be used to 


model W 


ij’ 


this would give 


W iJ ' - e ijk 3^ (2V T S kH> 


(2.15) 


where S 


kfc 


and 


V T are defined by Eqns. (2.12) and (2.14) > respectively. 


We shall refer to this as Model OJ-1. 


2.3.2 Model 03-2 


We note that the model given by Eqn. (2.15) involves computing the 


strain-rate tensor S 


k&’ 


which is an expensive process. It also uses 


the velocity field and hence requires storage space for the velocity fields 
even after the convective and stretching terms have been computed. Much 
computational saving could be obtained with a model that involves only the 
vorticity field; one such model is 


W 


ij 


if: <Vi> + s: <Vj> 


(2.16) 


where is defined by Eqn. (2.14). We shall refer to this as model 03-2. 

Both models U)-l and 03-2 can be shown to satisfy all three proper- 
ties mentioned previously (see Appendix A). Model 03-2 offers computa- 
tional as well as storage advantages over model U)-l and will be tested 
in Chapter 4 (along with model 03-1) , for the case of isotropic homoge- 
neous turbulence. 
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2.4 Filtering 

2.4,1 Sharp Cut-off (SCK) Filter 


Analytically, a filter that divides the large scales and the subgrid 
scales into two distinct regions in the Fourier sense would be convenient. 
Then f would contain all scales larger than a cut-off scale, and the 
subgrid scales (f') would contain all scales smaller than this cut-off 
scale. A one-dimensional version of such a filter is 

sin[k (x-x 1 )] 

G(x-x') = tt (2.17) 

7T(X-X ) 

and its Fourier transform is 

(0 if |k| > k 

H(k) = j c (2.18) 

( 1 otherwise 

We shall refer to this as the SCK (Sharp cut-off in k-space) filter. 

In inhomogeneous flows with turbulent rotational regions and irrota- 
tional regions, the two regions are separated by a sharp vorticity jump. 

In order to evaluate the ability of the SCK filter to smooth out jumps in 
the vorticity field, we apply it to a point vortex situated at the origin: 


w(x,y) = <5 (x) 6(y) (2.19) 


and 


w(x,y) 


sin[k c :xj sin[k c y] 
ttx iry 


( 2 . 20 ) 


W is plotted in Fig . 2.1. 

First we note that this filter creates oscillations and negative vor- 
ticity, which are undesirable from a physical point of view . Second , 
those oscillations decay slowly (they go as x ) , so the spreading into 
the irrotational region is excessive. 
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2.4.2 Gaussian (GS) Filter 


Another filter that has been used by previous workers (Kwak et al., 
1975) is the Gaussian spatial (GS) filter: 

G(x-x') “ ^ exp{-Y(x-x' ) 2 /A 2 } (2.21) 

where Y is a constant and A is the filter width. 

Applying this filter to a point vortex situated at 

get 

w(x,y) = j -4r exp{- (x 2 + y 2 )} 

11 tr b 

w is plotted in Fig. 2.2. 

We note that in this case we have created neither oscillations nor 
negative vorticity. By filtering the point vortex (Eqn. (2.19)), we have 
created another vortex with a Gaussian core of width A. 

We conclude that a Gaussian filter smoothes out jumps better than the 
sharp cut-off filter. Therefore, the GS filter was used in the cases in- 
vestigated in this work. 


the origin, we 

( 2 . 22 ) 


2.5 Computing Velocity Field from the Vorticity Field 

When the vorticity equation is used, the velocity becomes a diagnostic 
variable; i.e., the time variation of the velocity is not given explicitly 
by the equations but can be deduced once the vorticity is known. To do so, 
we shall define a vector potential ^ (see Lamb, 1932) such that: 


U i £ ijk 9x. ^k 
J 3 


(2.23) 


can be chosen to be solenoidal; i.e.. 


wM - 0 


( 2 . 24 ) 


Taking the curl of (2.23) and using (2.24), we get 


-,T ^ Ip. = “ 0). 

3x.dx. i i 

J 3 


(2.25) 
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Solving the Poisson equation (2.25) and using (2.23), we get the 
velocity field from the vorticity field. 

Note that the velocity field could have been deduced in another 
fashion by setting 


U). 


'ijk 3x. u k 
3 


(2.26) 


then taking the curl E p q^ °f Eqn. (2.26) to get: 

2 

__r_- ■ _ d ~ 

3x.3x. U i E ijk 3x. k 

3 3 3 


(2.27) 


and finally, solving the Poisson equation (2. .27), we get the velocity field, 
This approach involves differentiation of che vorticity field, followed by 
a double integration, whereas the first approach of (2.25) involves double 
integration followed by differentiation (2.23). Numerically, the first \ 
approach is usually more desirable; but in our case the two approaches are 
equivalent. Eqns. ( 2 . 2 3) — ( 2 . 25) will be used in this study. 

2.6 Summary 

Neglecting the molecular viscosity, the filtered dynamical equations 
in vorticity form become 


H . 3 

- A ~ — (u M . - u .0) . ) 
3t 3x. i x i l 

d 2 ip 

9x. 3x. W i 

3 3 


U i E ijk 3x. % 
3 


W. 


3x . ij 
3 


(2.28) 

(2.25) 

(2.23) 


where W . . is modeled as 


W ij ~ £ ijk 3x £ (2v T S kP 


(2.15) 


or 



or 


W. . 


3S7 ( Vi> + £: ( Vj> 

3 1 


( 2 . 16 ) 


where 


(C A) 2 (w.w.)^ 
v 11 


( 2 . 14 ) 


and 


S - — 

ij 2 


(^- u i + ^ u 3 ) 


( 2 . 12 ) 


with 


G(x-x') 


■ -• I 


T~ ex P \ “Y 
2^3 


(Xi-x |) 2 (x^x*) 2 (x 3 -x’) 2 


( 2 . 29 ) 


and 


A = (A 1 A 2 A 3 ) 


1/3 


It is in this form that the problem will be solved numerically. 


Chapter 3 
NUMERICAL METHOD 

Analytical solutions of the governing equations discussed in the 
previous chapter can be found for only very special cases, none repre- 
senting turbulence. Therefore, we propose using large computing machines 
to solve these equations for particular cases of interest. Numerical 
approximations of the governing equations require special care. In this 
chapter we discuss these approximations and present the methods we use 
to solve the difference approximation to the governing differential equa- 
tions . 

3.1 Notations 

A region of continuous space is divided into a uniform rectangular 
mesh; h. (i=l,2,3) represents the mesh width in the i fc direction. 

The mesh width need not be the same as the averaging width introduced in 
the previous chapters; we have used = 2tu and y = 6. For details 
on the effects of the filter width on the computational results see Kwak 
et al. (1975) . 

We then write the ^-component of the filtered flow quantity f- 

fch . 

at the n time step as 

4 n) = 1.2,3 (3.1) 

where (i , j , k) are the mesh point index for (x,y,z). 

We define the operator notation S/6£ to be the numerical approxi- 
mation to the continuous derivatives (9/3£) . 

3.2 Numerical Approximation 

Once space is discretized into mesh points, it remains to approximate 
the partial derivatives in terms of the values of the functions at those 
points. We have used two different approximation schemes: a fourth-order 
scheme and a pseudo-spectral method. 
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3.2.1 


Fourth-Order Scheme 


Using Taylor series expansions one can easily show that the approxi- 
mation to the partial derivatives, 

- t £- fu(i-2) - Mi-1) + 8u(i+l) - u(i+2)) (3.2) 

12h l 

is fourth-order accurate, i.e., the error in this approximation is of 
O(h^) . (For simplicity, the arguments j and k are not shown.) 

If periodic boundary conditions are to be used, u can be represen- 
ted by a discrete Fourier expansion (see next section) . 


u = ]£^u(k) e 1 — - 


(3.3) 


n 


where, for i =1,2,3, 


k i = 


n. 

1 


2tt * 

— — r— n . = wave number in the x . direction 

N ± h i x x 

N. ■ N, 

■ n -i _ 1 

~ , . ■ . , U, . . . , rj 4- 


N 


^ = number of mesh points in x^ direction 


u(k) is the discrete Fourier transform of u. Taking the discrete 
Fourier transform of (3.2), we get 


6u 

<5x. 


1 ( " 12h l k l n l2h l k l) g 

e - 8e + 8e - e > u 


12h 


{8 sin(h 1 k 1 ) - sin(2h 1 k 1 )} u 

Oil- 1 -L L ± 


where 


= ik|u 


k* = 
K 1 


gjp {8 sin(h^k^) - sinUh.^)} 


(3.4) 


is called the modified wave number. 


18 



Representation (3.4) allows us to evaluate the numerical approxima- 
tion (3.2) for the range of wave numbers up to rr/h^, the highest wave 
number that can be represented on a grid of size h . The Fourier trans- 
form of the exact derivative is ik.jU, so that, by comparing the modified 
wave number kj with k^, we see h6w well the approximation works (see 
Fig. 3.1). 

A more important consequence of representation (3.4) is that it 
allows us to integrate numerically in a manner consistent with our differ- 
ence approximation. In order to make this point clear, suppose we know 
the value, f, of the numerical approximation of the differential equation 



and we would like to find u, which when fourth-order finite differenced, 
gives us f exactly (to machine round-off). One way to do this is to 
write 

Au(i) « f(i) (3.6) 

0 8 -1 0 0 1 - 8 “ 

8 0 8 -1 0 0 1 

1 -8 0 8 -1 


10 0 

8 -i o ~ 8 ; °J 

for the case of periodic boundary conditions. This system of equations can 
then be solved in some standard way. 

Another way to handle this problem is by taking the discrete Fourier 
transform of . (3. 5)' to get 

ikju = f (3.7) 


where 


A = 


12h„ 


r 


Then, by solving for u, 


(3.8) 


multiplying (3.8) by e , and summing over all Jc, we obtain u. 

In this case only the one-dimensional transform is needed. This method, 
which is much more powerful than the previous one when integration in more 
than one direction is needed, will be used extensively for the solution of 
the Poisson equations (2.21). 

3.2.2 Pseudo-Spectral Method 

Periodic boundary conditions 

Suppose f(x,) is periodic in the x^ direction with period L (in 
the following we shall consider the one-dimensional case) and satisfies 
the "Dirichlet condition", i.e. , 

• f (x^) is defined at every point of the interval 0 x^ _< L, 

• f (x^) is everywhere single-valued, finite, and sectionally continuous, 

• f(x 1 ) is of "bounded variation", i.e., f(x^) does not have an infi- 

nite number of maxima and minima . 

It can be shown (Lanczos, 1956) that a function of this type can be 
expanded in a convergent Fourier series. 


f(x 1 ) = fOt-j.) e 

0^=-°° / 


■^^i x i 


(3.9) 


where 


T n i 


n» = — °°, • . . » 0 , 1 , . . . , 


f( k l) = t 


” ik l X l 

f(x 1 ) e dx x 


(3.10) 
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Since computers cannot handle infinite series, we have to truncate 

(3.9) , This is justifiable if f^k^) falls off rapidly for large | k.J ; 

this is the case of interest, since we filter the flow variables. Also, 
as mentioned before, we need to discretize in space. If mesh points 

are used in the x^ direction, the discrete analogs of Eqns. (3.9) and 

(3.10) become: 

N,/2-l 


where 


fUj) 



fOcp e 


lk l x l 


<3. XI) 


2v 

N,h, n l 


' 1“1 
*1 = jh l 


n. 


N N 

— 0 — - 1 
9 y •- • * > v * • • • 9 o x 


0, ... , N x - 1 


L/N, 


f(k x ) 


%-l 


' K T, f(x l> e 


- ik l x l 


3=0 


Fast algorithms (for = 2 n ; n = 1,2,...) have been developed 
(Fast Fourier Transform — FFT) by various workers (Cooley and Tukey, 1965 
Singleton, 1967) to evaluate the series (3.11) and (3.12) for the inverse- 
transform and the forward-transform, respectively. These will not be de- 
scribed in this work (we used a routine developed by Singleton, 1967). 

If we regard the expansion (3.11) as an interpolating formula, so that 
we treat 
tion, we obtain 


(3.12) 


x^ as a continuous variable, and differentiate the entire equa- 


tSf 

5x, 


“ iC f ^ k 1 ) ikj 

n. 


i^ C 1 X 1 


(3.13) 


The expansion (3.13) can be considered an approximation to the partial 
derivatives. Thus, to compute the partial derivatives of u, for the case 
of periodic boundary conditions, we proceed as follows: we find the dis- 

crete Fourier transform of the function in the direction in which the par- 

A 

tial derivative is needed, i.e. , we compute f(k^) from f (x^) ,■ ■' 


Multiplying f(k^) by ik-^ e and summing over all k^, we obtain 

6f/6x^. This is called the ’’pseudo-spectral” approach. This method has 
been analyzed by Lanczos (1956) and, with the development of techniques 
to compute the summations (3.11) rapidly, it has been proposed by Kreiss 
and Oliger (1973) as an approximation method and advocated by Orszag 
(1973) and Fox and Orszag (1973). 

For the range of wave numbers that can be captured with a given spac- 
ing and number of grid points and for periodic boundary conditions, the 
pseudo-spectral method yields extremely accurate values of the partial 
derivatives (see Fig. 3.1). 

The above method is limited to the case of periodic boundary condi- 
tions. However, the idea can be applied to other types of boundary con- 
ditions by using a set of orthogonal functions appropriate to the given 
boundary conditions. 


f = 0 boundary conditions 

If f(x^) is required to vanish at the boundary, i.e., f(x^) = 0 

for x^ = 0 and x^ = L, and is twice differentiable (a physically rea- 
sonable assumption), the Hilbert-Schmidt theory shows that its Fourier 
sine series 


f(x i> - £ v sin 


n^=o 1 
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L X 1 
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where 
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is absolutely and uniformly convergent. As in the previous section, we 
shall use the discrete analogs to (3.14) and (3.15), i.e. 


N l“ 1 


f(x 1 ) = ^ f (n^) sin 


n^=o 


n i 77 


(Nj-l)h^ "1 


K r 1 


f ' <n i ) " (i^T) ]T £<x i > sln 
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where 


* O f • • • f "l 

h l " VCNj-D 

x j^ = j « j * O, » o , - 1 

A s 

and f (n^) is the Fourier sine transform of f(x^). By using the 
FFT routine, a technique to compute the summation in Eqns. (3.16) and 
(3.17) can be rapidly developed. A detailed development of the Fast Dis- 
crete Sine Transform (FDST) is given in Appendix B. Generally, the FDST 
requires twice as much computation (for a given number of mesh points) as 
does the FFT. 


If we regard the expansion (3.16) as an interpolating formula, treat- 
ing x^ as a continuous variable, and differentiate, one obtains: 
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N-j-1 


<5x7 = £ fS(n l } k l cos 

1 n l=° 
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(3.18) 


where = n-jjr/(N^-l)h^ . In order to be able to use (3.18) as an approxi- 
mation formula for the partial derivatives, we need an FDST to find f (n^) 
we also need a Fast Discrete Cosine Transform (FDCT) . The discrete Fourier 
cosine series is defined in analogy to (3.16): 
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(3.20) 


where f'(n^) is the Fourier cosine transform of f(x^). 

Note that in (3.18) f S (0) * f 3 (N. -1) “0, making (3.18) exactly 

a s 

a discrete cosine transform of f (n^) . 

By using the FFT routine, a technique of computing the summations in 
Eqn. (3.19) and (3.20) can be rapidly developed. A detailed development 
of the Fast Discrete Cosine Transform (FDCT) is given in Appendix B. 

Thus, to compute the partial derivatives of a function which is 
zero at the boundary, we find its discrete sine transform f (n^) , mul- 
tiply it by n^Tr/ [(N^ - l)h^] and inverse transform using an FDCT routine. 
This method yields an extremely accurate approximation of the partial 
derivative when the function vanishes at the boundary, but its use is re- 
stricted to cases with a uniform mesh, 

9f/9x = 0 boundary conditions 

If f(x 1 ) is our function whose partial derivative Bf/Bx.^ vanishes 
at the boundary, i.e., Bf/3x 1 = 0 for x 1 = 0 and x 1 = L, then, by 
using arguments similar to those used before, it can be shown that its 
Fourier cosine series 

f(x 1 ) 

where 

f 

n. 


is uniformly and absolutely convergent. 

Equations (3.19) and (3.20) are the discrete equivalents of the above 
equations. If we regard expansion (3.19) as an interpolating formula treat- 
ing x^ as a continuous variable, then differentiate, we obtain: 

= M - 8S k * sin 


v 


(N^-l)h^ 


(3.23) 


OO 

- S 


n^=o 1 


cos 


n^TT 

IT X 1 


(3.21) 


/ f(x i 


) cos 


V 

L X 1 


dx. 


^c 

f ' 
n 


A C 

•f 

n. 


* "i * 0 


1 fi A c 

i\ • "i = ° 


(3.22) 


Obviously (3.23) satisfies the conditions 9f/9x^ * 0 at x^ * 0 
and x^ * L, and (3.23) is the discrete sine expansion of the partial 
derivative. 

Thus, to compute the partial derivative of a function for which 

/\C 

9f/9x^ =0 at the boundary, we find its discrete cosine transform f (n^) , 
multiply it by -k., and take the inverse transform using an FDST routine. 

The three methods described in this section will be used extensively 
as our approximation tools. 

3.3 Time Differencing 

To advance in time, a second-order Adams-Bashforth method was used. 

This method has been used by previous workers (Kwak et al., 1975; Shaanan 

et al., 1975), and use of a higher-order method was not felt necessary. 

If 3a)./3t = M., the Adams-Bashforth formula for a), at time-step 
i l i 

n + 1 is h 

. •. • • : / ....... • • • . • ■ 

^ + At || \ M^ n_1 ^ (3.24) 

In our case, 

3 _ z 8 

i 3x. ' ] i i 3 dxT ij 

Mote that this is a two-step explicit method. It is started with the 
Euler method : 

= 0)° + At M< 0) (3.25) 

3.4 Conservation Properties 

As was pointed out by Phillips (1959), numerical integration of the 
finite-difference analog of the Navier-Stokes equations may introduce non- 
linear instabilities if proper care is not taken. Arakawa (1966), working 
with the two-dimensional vorticity equation, showed that by properly con- 
serving vorticity, energy, and enstrophy (w^uk ) , these instabilities 
disappear. Lilly (1965), working with the primitive variables, developed 
a spatial-differencing scheme that conserves momentum and energy. By 


conservation we mean that, in the absence of external forces and viscous 
dissipation, the only way that the momentum and kinetic energy in a con- 
trol volume can change is by flow through the surface. This property must 
be retained by the numerical approximation. In the simple case of periodic 
boundary conditions, we have 
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dt j 
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o 

(i.e., momentum conservation) 

(3.26) 


) 2 u iV iv ■ 0 

(i.e., energy conservation) 

(3.27) 

usually easy to devise a 

numerical approximation to the 

dynamical 


equations in primitive form that conserves momentum, i.e., summation over 
the flow volume of the approximate equations would give the discrete 
equivalent of Eqn. (3.26). However, the difficulties arise when trying 
to show energy conservation, since in general the identity 


u . 


i 3x. u i 
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3 1 

3x . 2 u i u i 
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(3.28) 


does not hold in finite-difference form. 

Writing the equations of motion in the following form (Tennekes and 
Lumley, 1972) : 
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and integrating over the flow volume, we get 
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f u dv+ f u . (t± . dv . f » (i + i u ,,.) dv 

J6- X 3 3x ± y 2 3 3 ! 

For periodic, boundary conditions, integration by parts yields: 

f u. u. dv = - f u. u. dv = 0 (using (3.30)) 
JU- 3 3x. i x ox. 3 
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and 


+ * u J u j) 


dv = 0 


Hence Eqn. (3.31) reduces to Eqn. (3.26) and we have momentum conservation. 
Now, multiplying Eqn. (3.29) by u^, we get: 
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9t 2 u i u i 
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where the convective terms sum to zero by symmetry. 

Integrating (3.32) over the entire domain yields: 
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at 
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For periodic boundary conditions, integration by parts yields: 

f a /p , i \ y fi? V i \ 9 , 

/ u -t r + ^ u. iJ. I dv = ~ / [■ T + “T u. u. ) -5 — u . dv 

^5- i 9x i \ P 2 3 3 / y \ P 2 3 J / 9x i 1 


dv = 

= 0 (using (3.30)) 


Hence Eqn. (3.33) reduces to Eqn. (3.27) and we have energy conservation. 

We notice that, with the equations written in the form Eqn. (3.29), 
we did not need the identity (3.28) to show energy conservation from the 
dynamical equations in primitive form. The conservation properties were 
obtained by making use of only integration by parts and the continuity 
equation. ■■■■ 

Consider the numerical approximation of Eqns. (3.29) and (3.30): 
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where we are using 6/6x^ to denote the numerical approximations to the 
partial derivatives 9/9x^, and the same approximations are used in both 
equations (3.34, 3.35) for any given independent variable. In order to 
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have long-term Integration stability, Eqn. (3.34) should numerically con- 
serve momentum and energy. 

If we follow the steps used in deriving the conservation properties 
from Eqns. (3.29) and (3.30), we realize that the conservation properties 
will follow if we can establish numerical summation by parts. Consider 
the oiK^-dimensional case, where we have, for periodic boundary conditions, 


/■ 


(x) f(x) dx 




f(x) u(x) dx 


The numerical analog of the above equation is: 
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Expanding u(j) in Fourier series, we get: 

u(j ) = 2^ u(n) exp(2irijn/N) ; j = 0,1, . . . ,N-1 
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where the u(n) are given by the inverse transform: 
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(3.37) 


where k'(n) is the modified wave number. The modified wave numbers for 
the numerical methods we are using are 

• ik' = ik for pseudo-spectral, (3.38) 


• ik' « 1 gh sin ( ktl ) “ sin(2kh)] (fourth-order approximation. 

Substituting Eqn. (3.37) into the left-hand side of Eqn. (3.36) yields 

N=1 r -j N4 N/2^1 

«U) • 6 1 ^ ^ ^ 
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f(j') 
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Now, changing the summation index in the last sum from n to -n, we 
see that this expression will agree with the right-hand side of Eqn. (3.36) 
provided that: 


k ? (n) = - k'(-n) 


(-!) - o 


(3.39) 

(3.40) 


Condition (3.39) is satisfied by all the methods under consideration, 
and k'(-N/2) *0 is true for the finite-difference method. The pseudo- 
spectral method cannot differentiate between f = exp(ijir) and f = 
exp(-ijTT), and, due to this confusion at n = -N/2, k'(-N/2) is set 
equal to zero for the pseudo-spectral method. Hence, summation by parts 
is obtained when (3.39) and (3.40) hold. Summing Equation (3.34) over 
all mesh points, using the generalization of (3.36) to three dimensions and 
using Eqn. (3.35) yields the numerical equivalent of (3.26). Multiply- 
ing Eqn. (3.34) by u^, the nonlinear term in the left-hand side of (3.34) 
will sum to zero by symmetry; then, using as before the three-dimensional 
generalization of (3.36) and (3.35), summing over all mesh points will 
yield the numerical equivalent of Eqn. (3.27) . 


3.5 Differenced Vorticity Equations 

In order to insure that the numerical approximation to the vorticity 
equations are equivalent to the numerical primitive equations, we must 
take the numerical curl of Eqn. (3.34). Before doing so, we note that, 
numerically. 


V • V X v 
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(3.41) 
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where V and S are any vector or scale, respectively, the above ex- 
pressions are identically zero, if for each direction the same approxima- 
tion is used for all operators. 
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The numerical curl of (3.34) is 




(3.43) 


Equation (3.43) conserves vorticity, i.e., summing it over all space 
the total vorticity in any control volume (subject to periodic boundary 
conditions) does not change with time. Hence in the form (3.34), the 
primitive equations also conserve vorticity. 

The numerical divergence of (3.43) is 


3 fi 


at 6x. i 

r 


u> . = 0 


Therefore, an W field solenoidal at time t will remain solenoidal at 
time t + At. 

3.6 Poisson Equation 

Having the vorticity field ok at time step n, we have to find the 
velocity field in order to be able to advance in time. To do so, we shall 
define a vector potential (also called the vector stream function) i/^, 
such that 




u i ■■■• £ ijk fix: r k 


can be chosen to be solenoidal; i.e.. 
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fix. y i 
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Taking the curl of Eqn. (3.44) and using Eqn. (3.45), we get 


4 >, = - a) , 


fix. fix, ^i 
J J 


(3.46) 


The Poisson equations (3.46) will be integrated using the approach 
introduced in Section 3.2. For the case of periodic boundary condition, 
the discrete Fourier transform of Eqn. (3.46) is 


ktkt \p. 
3 3 


= - W, 


(3.47) 
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where k! is the modified wave vector introduced in Section 3.2, 

A 1 

for we have 


Solving 


(3.48) 


and by inverting the transform we obtain the stream vector consistent with 
our numerical differencing. It satisfies two conditions. First, the vel- 
ocity field obtained using (3.44) will be solenoidal. We have in Fourier 
space: 


k! u. 
l i 


E yk k i k j ♦k = 0 


(3.49) 


Second, taking the curl of (3.44), we have in Fourier space: 


ok = e . .. ik! u. = - e. E. k!k' 

i ijk j k ijk kpq 3 P <1 


(3.50) 


= — k| k! lb . + k!k! \p. 

T i T l X 3 r 1 


Since k! liJ. = 0, (3.50) is exactly the Poisson equation (3.47). 

j J 


з. 7 Numerical Filtering 

Examination of Eqn. (3.24) reveals that the on ly numerical problem 

left is the numerical evaluation of the u.td . - u .to. term. Since u .to. 

4 -.**• .... 4 4 

и. w. can be computed easily, the problem is that of numerical filtering. 

Filtering is the evaluation of a convolution integral 


u.to. G(x-x') dx* 
3 1 


(3.51) 


If this integral is evaluated using conventional integration routines, the 
computation cost is prohibitive. Previous workers (Leonard, 1973; Kw.ak et 
al., 1975; Shaanan et al., 1975) argued that the filtered terms u^(x’) 
and u (x») are smooth, and they expanded those terms in a Taylor series 
about x. Using a Gaussian for G(x) , they obtained: 


A 2 r-,2 


u j U) i = uyo ± + ^ V (u .to ± ) + 0(A ) 
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and the 0(A ) term was called the Leonard term. The above approxima- 
tion will require the use of a fourth-order, finite-differencing method 
(Kwak et al, 1975) or a modified second-order method (Shaanan et al., 1975) 
that yields the Leonard term as its truncation error. However, when 
higher-order methods are used the expansion (3.52) needs to be extended to 
higher orders, and the computational expense becomes prohibitive. When 
periodic boundary conditions are used, we can take the Fourier transform 
of Eqn. (3.51) to get: 


u.O). = (u.io.) G 

J i J i 


(3.53) 


Thus, given u^ and 0)^ , one can compute the term (u.U)..), multiply 
it by G, then simply invert the transform to obtain u^UK. 

When u.O). vanishes at the boundaries, i.e. , u.O). =0 at x = 0 

1 i 1 i 

and x = L, we can expand it in a Fourier sine series. Taking the one- 
dimensional case, for simplicity, we set 

AS 


u j“i • £ sin (t x ) 
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n=o 

Substituting (3.54) in (3.51), we get 
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u.o) ± = J 2^ (u.^) sin (x— x')) G(x’) dx* 


Since the series (3.54) is absolutely and uniformly convergent, we 
can take the summation outside the integration to obtain 
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UjOJi = u j w ^| jf sin(~- x) cos - x’) G(x') dx' 
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If G(x) is an even function, which is the case when (2.21) is used, 
the second term in the bracket vanishes and one obtains 
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where 


= J' G(x’) cos (^~ x’j dx' 


is the Fourier cosine transform of the Gaussian filter. 

What Eqn. (3.55) tells us is that, for the case in which u.U)^ = 0 

for x = 0 and x = L, u.U). can be computed by the following procedure: 

3 1 

we first compute the Fourier sine transform (u.U)^) of u U)^, and then 

A ^.3 J 

multiply it by the Fourier co sine transform G of the filter, to obtain 

the Fourier sine transform (u_.0)^) of u^U)^. Finally, inverting the sine 

transform, we obtain u.U).. 

3 i 3 

Similarly, it can be shown that, for the case in which u^uh = 0 
at x = 0 and x = L, we have 


u i W i “ S G cos (t" x ) 

J n=o J 


(3.56) 


or 


(u.co.) 
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(u.U).) G 
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By the use of the FFT, FDST, and FDCT, ’'exact" filtering can be ob- 
tained for all boundary conditions of interest with acceptable computa- 
tional speed. 

An important property required of a filter is that the filtered value 

of a constant must be the same constant. Numerically, it is desirable to 

preserve this property, which is equivalent to requiring the integral of 

/\ 

the filter function be unity or G(0) -1. The exact continuous Fourier 
transform of (2.21) ,is-. 


G(k) = exp^- ~ k 2 ) 


(3.57) 


When G(k) is discretized, we get 

S D (k) - ex p(~ (I 1 n) 2 ) . n = 0.±1,±2, 


(3.58) 


Hence G^(0) = 1. 


Another property required of a filter is that it smooth out jumps 
(see Section 2.4) without introducing oscillations. We have modeled the 
situation with a top-hat function: 


f(x) = 

Analytically, we have 


t 1 

0 


x x < x < x 2 
otherwise 


(3.59) 


f(x) > j (erf (x^ - x) - erf(x 2 ~x)) 


which is a smooth function with no oscillations. 

A 

When (3.59) is discretized and filtered numerically using G^(k) , 
the top-hat function, Eqn. (3.59), is smoothed out (see Fig. 3.2). How.- 
ever, small oscillations are introduced. This is due to the fact that the 
discrete inverse transform of (3.58) is not smooth. For this reason we 
have used a discrete Gaussian in x-space. 


where 

where 


; D (x> - iexp^-Y^ 
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x = hn 


■ ? exp (- 


(3.61) 


as our filter function. The oscillations in the x-space (see Fig. 3.2) 
do not appear when this filter is used. 


3.8 Summary 

The dynamical equations in vorticity form will be solved as follows: 
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Chapter 4 

DECAY OF ISOTROPIC TURBULENCE 

4.1 Background 

In order to assess the feasibility of using the vorticity equations 
as the governing equations for turbulent flows , we applied the computa- 
tional methods described in Chapter 3 to the simplest problem in turbu- 
lence, namely, the decay of homogeneous isotropic turbulence. This flow 
was also used to determine the value of the subgrid scale model constant 
for use in subsequent calculations of other flows. 

The grid turbulence experiment of Comte-Bellot and Corrsin (1971) was 
used as the "target" for our numerical predictions. When viewed in a co- 
ordinate frame moving with the mean velocity, this experiment approximates 
homogeneous isotropic turbulence. 

This study was presented in an earlier report (P. Moin et al., 1978) 
and is rediscussed in this work to support the argument that model to- 2 used 
in conjunction with the vorticity equations produces similar results to 
those obtained using the more commonly used model to-1 . The contributions 
of Hr. P. Moin are gratefully acknowledged. 

4.2 Initial Conditions 

We started with an initial field that is divergence-free and has a 
spectrum obtained by filtering the experimental spectrum at the non- 
dimensional experimental time T = U Q t/M - 42. U Q =10 cm/sec is the 
experimental free-s tr earn air speed, M = 5.08 cm is the size of the ex- 
perimental turbulence-generating grid, and t is the real time in seconds. 
The initial field was otherwise random. The generation of such a field is 
discussed in detail by Kwak et al. (1975) and will be briefly outlined 
herein. The filtered field is generated in k- space by setting: 

“i® = (T) 3 (^) a ( aA i + lbB i> <«■» 

where 
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¥(k) = filtered experimental energy spectrum at time T » U Q t/M = 42, 

a = cos(0) 
b = sin(0) 

where 0 is a random angle, and are unit vectors picked such 

that A^k^ = B^k^ = 0, otherwise random. 

To insure that (4.1) is the Fourier transform of a real field, we 
must have 

u ± (k) = u*(-k) ' (4.2) 

/S 

where * indicates complex conjugate. Now, by inverse transforming u ± , 

u . is ob tained . 
x 

Using the above initial field, we shall use the methods of Chapters 
2 and 3 to predict the spectrum at T = 98. The predicted spectrum will 
be compared with the filtered experimental spectrum at T = 98. 

4.3 Selection of C„ 

The model constant was obtained by matching the computational rate 
of filtered energy decay to that of the experiment (Fig. 4.1). The values 
of the constants obtained using different numerical schemes and different 
models were in most cases within ten percent of each other (C^ = 0.2 + 0.02 
see Table 4.1). 

4.4 Results 

Under the assumption that the computational box size is large compared 
to the scale of the energy-containing motions, we can use periodic boundary 
conditions in all three directions. A uniform cubic mesh system was used 
with N, the mesh number in each direction, and h, the mesh spacing, 
chosen such that the computation captures as much of the turbulence energy 
as possible (Kwak et al. , 1975), We used the sets 

■■V' 1 16 h = 1.5 cm , t = 6.25 x 10 sec 

and 

' ■ ■■ -3 

N = 32 , h = 1.0 cm , t = 6.25 x 10 sec ; ^ 
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When periodic boundary conditions are used, it was shown in Chapter 3 
that the pseudo-spectral method is more powerful than any finite-difference 
method. However, when the periodic pseudo-spectral methods cannot be used, 
we may have to use finite-difference methods. Since one of our objectives 
is to determine the model constant for the vorticity equations, both the 
fourth-order finite differencing and the pseudo-spectral methods were ap- 
plied to the case of isotropic homogeneous turbulence. 

4.4.1 Fourth-Order Finite Differences 

Figure 4.2 shows the energy spectrum obtained by fourth-order finite- 

differencing the vorticity equation, using model to - 1 (Eqn. (2.1 )) for 

3 

the subgrid-scale model, on a 16 mesh. Our results compare well with 
the experimental results up to wave number 2.5, after which the inaccu- 
racy of fourth-order differencing begins to show. Fourth-order differenc- 
ing the primitive equations (Kwak et al, , 1975; Moin et al., 1978) produced 
good agreement with the filtered experimental results using the primitive 
variable version of this model. This shows that the vorticity approach is 
equivalent to the primitive variable method. Thus the use of the vorticity 
equations is definitely feasible in turbulent flow computations. 

4.4.2 Pseudo-Spectral Method 

Figures 4. 3-4. 6 show the energy spectra obtained using the pseudo- 

3 

spectral method, with 16 mesh. Fig. 4.3 shows the results obtained using 
model to- 1 (Eqn. (2.15)). We note that for k > 1 the computed results 
are considerably lower than the experimental values. This indicates that 
the subgrid-scale model is draining too much energy from the small struc- 
tures, and, since our total energy is equal to that of the filtered experi- 
mental value, too little energy is taken out from the large structures. In 
this case, we used the pseudo-spectral approximation to calculate the sub- 
grid scale terms as well as the other terms. 

Figure 4.4 'shows the energy spectrum obtained using second-order cen- 
tral differencing to approximate the derivatives appearing in the subgrid- 
scale model (see Section 3.. 8) : 
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61 f(i+l) - f (1-1) 

6x 2h 

We note a considerable improvement in the spectrum, except for a small 
accumulation of energy at the extreme (high wave number) end of the spec- 
trum, which was present to a lesser extent in Fig. 4.2, / 

3 ! 

Figures 4.5 and 4.6 are the results from a 16 computation using 
the pseudo-spectral method and model 0) - 2 (Eqn. (2.16)) for the subgrid- 
scale model. We note the same behavior in Fig. 4.5 as in Fig, 4.3; the 
computed spectrum falls below the experimental spectrum, indicating that 
using the pseudo-spectral method to compute the spatial derivatives in 
the subgrid-scale model damps too much energy in the wave number range 
k 1, Using second-order finite differencing to compute the partial 
derivatives in the model to- 2 (Eqn. (2.16)) , we obtain a significant 
improvement in the computed spectrum (Fig. 4.6). These results are similar 
to the results obtained using model to-1, indicating that the two models 
are equally good. 

3 

Figure 4.7 shows the energy spectrum obtained from a 32 pseudo- 

spectral calculation, using second-order finite differencing to compute 

the partial derivatives in model to- 2. The results are similar to those 
3 

of the 16 computation. 

It can be concluded from these results that the vorticity equations 
provide a satisfactory basis for the simulation of homogeneous isotropic 
turbulence. Both models to-1 and to- 2 produce similar results. Model 
to- 2, given by Eqn. (2.16), will be used in the following computations, 
due to the computational advantages it offers over model to-1 (see Sec- 

3 

tion 2.3). Finally, a relatively coarse 16 mesh is sufficient to cap- 
ture interesting features of the homogeneous isotropic turbulence, and no 
significant improvement in the energy spectrum was obtained by using a 
32 ■ mesh system. , 

4,5 Computational Details 

The calculations described above were executed on the CDC-7600 at 
NASA-Ames Research Center, using programs written in Fortran. The total 
storage requirements (octal) were as follows : 


3 

16 Calculation 



Large Core Memory: 

Four th-order 

310,360 



Pseudo-spectral 

230,000 


Small Core Memory: 

Fourth-order 

104,465 



Pseudo-spec tral 

61,334 

32 3 

Calculation 




Large Core Memory: 

Pseudo-spectral 

1,110,000 


Small Core Memory: 


126,605 

The computing time per 

computational time 

step was approximately as fol- 

lows 

• 



16 3 

Calculation 




Fourth order - 2 

.5 sec CPU time 



Pseudo-spectral = 

4.0 sec CPU time 


32 3 

Calculation 




Pseudo-spectral = 

34 sec CPU time. 
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Chapter 5 

MIXING LAYER: TWO-DIMENSIONAL COMPUTATION 

5.1 Preview 

It Is well documented (Wlnant and Browand (W&B), 1974; Brown and 
Roshko (B&R), 1974; Konrad, 1976; Dimotakis and Brown (D&B) , 1976) that 
in some cases the spatially developing mixing layer contains coherent 
structures (in the terminology of B&R) or discrete vortices (in the termi- 
nology of W&B), In these experiments, the mixing layer grows via the 
interaction of neighboring vortex-like structures that rotate around and 
combine with each other to form a similar but larger structure (see Fig. 
5.1) . This mechanism is called vortex pairing. In this chapter we study 
the vortex-pairing mechanism by perturbing an infinite array of vortices. 
The effect of the initial perturbation on the roll-up is discussed. All 
cases treated in this chapter are completely two-dimensional; three- 
dimensional cases are discussed In the next chapter. 


5.2 Some Experimental Results 

The mixing layer is generated in a laboratory by bringing together 
two streams of fluid of different streamwise velocity (see Fig. 5.2). 

The measured mean velocity profiles, at different streamwise positions, 
are self -similar and can be fitted by an error- function (Spencer and Jones 
(S&J), 1971): 


* * 2 ~~ fe ~ erf (a(u-n Q ) 


where - 

r “ u 2^ u l’ 

u^ = velocity of the high-speed side, 
u„ = velocity of the low-speed side, 
r| = z/(x-x q ) 
a = spread parameter 
z - cross-flow coordinate, 


(5.1) 




X 


x = streamwise coordinate, and 
x Q = virtual origin of the layer. 


Rearranging (3.1) and normalizing the velocity on Au = 


U 1 ” u 2* we 


‘-ir - °- 5 erf(o(n-n o » 


(5.2) 


where U = + /2 is the mean velocity. The spread parameter O is 

a function of r, and the spread data can be fitted by the expression: 


(5.3) 


where CT q is the spread parameter for r = 0. S&J report a Q = H for 

other values of a : see Table 1.1. 

o 

Defining the momentum thickness, 0, to be 


1 

= J I (U •“ Un) (U-| ~ 

(Au)^y-oo 


u) dz 


(5.4) 


and substituting (5,2) in (5.4) , we get 


o Au 


a /2ir a 2 /2tt 

o 


(5.5) 


e m 


(5.6) 


Since O is constant, Eqn. (5.5) shows that the momentum thickness grows 
linearly with x. 

Substituting (5.6) in (5.2), we get: 


= 0.5 erf 


/ (z ~ z c>\ 

\ e 72 ? / 


(5.7) 


Due to computer limitations, one cannot set up a uniform grid that 
covers the length of the experimental set-up (1.8 m for the W&B case) and 


at the same time resolves the large-eddy scale (^>1—4 cm). We propose to 
use a uniform grid that moves with the mean speed U. The size of the 
computational domain is chosen so that two vortices are captured in the 
initial field; i.e. , we can imagine that we are following the fluid in 
the dashed box in Fig. 5.1 as it moves downstream. 

In our frame the layer will develop in time rather than space. We 
shall in fact be studying a portion of a time-developing mixing layer. 

This layer can be thought of as being created by having two infinite coun- 
termoving streams of velocity ± Au/2 brought in contact suddenly at T = 0 
For this flow, the mean quantities will be horizontal planar- averaged quan- 
tities; for example, the mean velocity profile will be defined as 


< u > 


xy 


//“ 


(x,y,z,t) dx dy 


(5.8) 


The momentum thickness, defined as 

0(t) 




(5.9) 


will be a function of time instead of space. According to the Taylor 
hypothesis, the state of the flow at the experimental streamwise distance 
x is the same as that of the computed layer at the computational time 
variable t. The variables x and t are related by the expression: 


x = Ut (5.10) 


Substituting (5.10) in (5.5), we get an expression for the expected 
momentum thickness of the t ime-developing layer: 


9(t) = 


t - t 


a 2 
o 


Au 


(5.11) 


Equation (5.11) shows that 0(t) should grow linearly with time, with 


d9 

Audt 


a 2 /2 tt 
o 


(5.12) 
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5.3 Boundary Conditions 


The coordinate system used is shown in Fig. 5.3, where the x-direction 
is the strearawise direction, the y-direction is the spanwise direction, and 
the z-direction is the cross-flow direction. We shall use periodic bound- 
ary conditions in the x- and y-directions; this is allowed if the size of 
the computational box is sufficiently greater than the integral scale in 
a given direction. At a large enough z location the flow is essentially 
horizontal and uniform. We can use no stress boundary conditions in the 
z-direction (i.e., 3u/3z = 3v/3z = w = 0 at z = 0 and z = L) if the 

boundaries of our box in this direction are sufficiently far from the cen- 
ter of the layer. This will allow us to expand the velocity fields as fol- 
lows: 


u = SS'2 u ( k 1# k 2 , n ) e 


n k 2 k 1 


n k 2 k x 


n k 2 k^ 


i(k^x + k 2 y) 

cos| 

(#} 

(5.13) 



V 3 / 


i(k^x + k 2 y) 

cosj 

Imrz \ 

1^7 

(5.14) 

iC^x + ^y) 

sinj 

(mrz\ 

V l 3 I 

(5.15) 


and the vorticity fields as follows: 


, 7 , (^(k-^k^.n) e 

n k 2 k^ 


w 2 “ V . / > / : -.v : (^(k^Vk^in) e 
n k 2 k x 


( 0 , 



11 k 2 k i 


a) 3 (k 1 ,k 2 ,n) e 


i(k 1 x + k 2 y) 

sinj 

1 n?r z \ 

c m 

(5.16) 

i(k^x + k 2 y) 

sinj 

i nTTz V 
\ L 3 / 

(5.17) 

i(k^x + k 2 y) 

cosj 

1 m\z \ 
/ 

(5.18) 


The pseudo-spectral method will be used to approximate the partial 
derivatives. The numerical technique was discussed in Chapter :3> 7 7 
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5.4 Initial Conditions 


We want to prescribe an initial profile that corresponds to a pair 
of vortices. It has been shown in Chapter 2 that filtering a line vortex 
produces a vortex with a Gaussian distribution of vorticity in the core. 

We shall use this fact to generate our initial conditions. 

The initial conditions are generated by starting with two line vor- 
tices in the spanwise direction at (x = x^, z = L 3 / 2 ) and (x = X 2 » z = 1 , 3 / 2 ) 
(see Fig. 5.4), and filtering in the x-z plane with the relatively wide 
Gaussian filter: 


G(x, z) = 


__i_ (A A\ 

A 1 A 3 \ 6 h^ 6 h 3 / 


(5.19) 


where h. is the mesh size in the i-th direction (i = 1,3) and A, 


(i = 1 , 2 ) 
field: 


is defined by Eqn. (3.61). This will produce the vorticity 


u>- 


» C 


1 A 1 A 3 


exp 





/ (z-L 3 /2) Z \ 

e T -*r) 


0 < x < L, 


0 < z < L, 


(5.20) 


w. 


= w. 


= 0 


£a) 2 ^ x 1 ,z ^ = w 2 ^ x n ^i» z ) n = 11 , ± 2 , .. . (periodicity) 


where C^ is an arbitrary constant that adjusts the strength of the vor- 
tices. Note that these vortices can be elliptical; they are 113/113 times 
as long in the streamwise direction as in the cross flow direction. 

Equations (5.20) correspond to a perturbed infinite array of vortices 
with a perturbation parameter 3 equal to: 



(5.21) 


3=0 corresponds to a uniform (unperturbed) vortex array, and we need 
to deal only with the case 3 > 0 . 
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Figures 5.6a-f show constant vorticity contours for $ = 6/16, 5/16, 
4/16, 3/16, 2/16, and 1/16. Note that for large 3 the vorticity con- 
tours look like those for a single distorted vortex. 

5.5 Mesh Size Selection 

We have shown in Chapter 4 that a 16x16 x16 mesh system can re- 
solve isotropic homogeneous turbulence with acceptable accuracy. For the 
cases considered in this chapter there are no variations in the spanwise 
direction. We dropped the number of meshes in the spanwise direction to 
^ = 4, the minimum number of meshes that our three-dimensional code was 
designed to handle. In the cross-flow direction the mesh number was in- 
creased to = 33 in order to allow the layer to grow in this direction. 
This gives a total number of mesh points of x x = 16 x 4 x 33 

= 2112 . 

The spanwise vorticity is defined by 

-k “ - m 7 (5 - 22) 

Averaging (5.22) over x-y planes and using periodic boundary conditions, 
we get 


< u>„ > 

2 xy 


-j- < u > 
dz xy 


(5.23) 


If we substitute in (5.23) the vorticity distribution given by Eqn. (5.20), 


- 7 - < u > = C 

dz xy 


2 / < 
1 exp ^- “ 


z-L 3 /2)' 


(5.24) 


This ordinary differential equation can be solved together with the bound- 
ary condition: 


< u > = 0 at z = Lo/2 (5.25) 

xy 3 

The solution is obtained by simple integration: 




■ .fee A.— fflSR i i’*» ■ ■ 


< u > 


i eJrVi) 

1 V6 h_ / 


xy L 


Non- dimensional! zing the velocity Au, we get: 


< u > = 0.5 erf 

xy 


/ ZL 3 /2 \ 

VS h, / 


Equating Eqns. (5.26) and (5.27) and solving for C^, we get: 


(5.26) 


(5.27) 


C 1 = 0,5 ^ 


(5.28) 


The length scales are non-dimensionalized on the momentum thickness. 
The mesh size was chosen such that the initial momentum thickness is equal 
to unity. Substituting (5.27) in Eqn. (5.9), we get: 


), = — h- = 1 

in /2? 3 


and solving for h^» we obtain 


h 3 


= 1.023 


(5.29) 


The mesh size in the streamwise direction was set equal to: 


h l = 


1.364 


(5.30) 


The non-dimensional time step was picked up to be equal to: 


AT ■= = 0.0799 

in 

which yields a Courant number such that : 


(5.31) 


U — < : 0.03 

C : . 00 h- “* 

which is well within the stability criterion and assures that the error 
caused by the time advancement will be acceptably small. 

The mesh size in the spanwise direction is irrelevant for the cases 
considered in this chapter. We have set h 2 = h^. 


5.6 Selection of 3 


We have shown in Section 5.2 that, to accord with the experimental 
observations, the momentum thickness 0(t) must grow linearly with time; 
and, using = 11 (S&J), we expect: 


d9 

Audt 


- » 0.018 
a 2 /2? 


(5.32) 


o 

We have run a series of calculations for different values of 3. 

Fig. 5.5 shows the momentum thickness plotted vs. T for the 

cases run. For the highly perturbed cases, 8 > 4/16, the momentum thick- 
ness 0(t) does not grow linearly in time. However, for 0 - 3/16, 2/16, 
and 1/16, 0(t) does grow linearly in time, with d0/Audt = 0.020, 0.015, 

and 0.009, respectively. 

Figures 5.6 and 5.7 show constant vorticity (contour) plots for the 
various cases at times T = 0 and T = 16.78, respectively. Figs. 5.6a-c 
and 5.7a-c show that for large 0 we have essentially one elliptical vor- 
tex which grows "fatter” in time, to become more or less circular at time 
T = 16.78. Figs. 5.6d-f show that for small 0, we have initially two 
distinct vortices; these vortices draw closer and rotate around each 
other (Figs. 5. 7d-f) . For the case 0 = 3/16, the two vortices merge to 
form one vortex at time T = 16.78 (Fig. 5.7d) . 

The above observations indicate that case 0 = 3/16 gives results 
comparable to the experimental observations. The spread parameter 
obtained for 0 = 3/16 is equal to 


o - 
o 


d0 

Audt 


2 >/2tt 


= 9.97 


which is within 10% of the experimental results of S&J. 


5.7 Mean Velocity Profiles 

The mean velocity profile < u > xy defined by Eqn. (5.8) is a func- 
tion of z and T. Fig. 5.8 shows 2< u > x ^/Au plotted vs. z/0 at 
AT = 2.4 Intervals, for 0 - 3/16. The profiles collapse into one, indi- 
cating self-similarity of the mean velocity profiles. 
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Self-similarity is also observed in the experimental data. Thus, as 
far as the mean profile is concerned, the data can be fit by pairing vor- 
tices with 3 = 3/16. 


5.8 Mean Turbulent Intensity Profiles 

In our computational box, the non-dimensional mean turbulence inten- 
sity is defined as 

2 

— 9 — = - — - < (.u — < u > ) ^ + (v - < v > ) ^ + (w - < w > ) ^ > 

2 xy xy xy xy 


(5.33) 


2(Au) 2(Au) ' 

where < > are planar averages defined by Eqn. (5.8). 

xy ' .. .. 

Figure 5.9 shows the mean turbulence intensity plotted vs. z/0, for 
the case 3 - 3/16, at At = 2.4 intervals. We note that the turbulence 
intensity decays slightly at the early stages of the pairing and then 


reaches a self-similar situation 
:ed 

i~2 


2 2 

Compared with the experimental results, our peak intensity q /2(Au) 


max 


= 2.06 x 10 is substantially lower than the experimental value reported 
-2 

by S&J (3.5 x 10 ). The low value of the maximum turbulence intensity is 
due to the fact that we did not take into account the subgrid scale contri- 
butions, and that our field is strictly two-dimensional, whereas in reality 
spanwise fluctuations are present in the experiment of S&J. 


5.9 Summary 

It is interesting to note that vortex pairing is capable of producing 
self-similar mean velocity and turbulence intensity profiles, and a linear 
growth of the momentum thickness that compare with experimental results 
(for 3 =3/16). We note that, due to periodic boundary conditions, once 
the vortices have paired we get a uniform vortex array (3 =0) and the 
pairing and layer growth stop. If we want the pairing to continue, We 
would have to perturb the array by displacing the vortices in the stream- 
wise direction. We have not done this because in the actual flow succes- 
sive pairings are not clearly separated and are random. 


A uniform array of vortices can be perturbed in several different ways; 
for example, by adding a cosine distribution of vorticity to a uniform array, 
we can enhance the pairing (see Appendix C) and get results similar to the 



results presented in this chapter. One could also make the vortices of 
different strengths or use any combinations of these perturbations. 

The perturbation 3 = 3/16 (vs. 3*0 for the unperturbed layer) 
needed to achieve the observed experimental growth rate of the shear layer 
may, at first, seem excessive. In the experiments, the downstream vor- 
tices exert a significant influence on those in the initial portions of 
the layer (D&B) ; recall that the influence of a distant vortical struc- 
ture on a given point decreases inversely with distance. The cumu- 
lative effect of the downstream vortices can be considerable, and, since 
they tend to be highly turbulent, they may strongly perturb the vortices in 
the initial section of the mixing layer. Therefore, the value 3 = 3/16 
may in fact be quite reasonable. 


Chapter 6 

MIXING LAYER: THREE-DIMENSIONAL COMPUTATIONS 

6.1 Preview 

In Chapter 5 we started with a two-dimensional initial field, and the 
numerical simulation of the governing equations stayed two-dimensional. 
However, actual flows are rarely two-dimensional, and truly turbulent 
flows are always three-dimensional. (Two-dimensional turbulence is ap- 
proximated by certain atmospheric structures and in highly stratified 
fluids.) In this chapter we evaluate the importance of large structures 
in the development of the mixing layer, which is two-dimensional in the 
conventional mean sense but contains the three-dimensional structures. 

6.2 Boundary Conditions and Mesh-Size Selection 

The boundary conditions and coordinate system of Chapter 5 will be 
used. Periodic boundary conditions will be used in the streamwise (x^) 
and spanwise (X 2 ) directions, and no-stress boundary conditions in the 
cross-flow (x^) direction. 

The number of meshes used for the cases discussed in this chapter is 
16 x 16 x 33 = 8448. The mesh sizes and time step are the same as in the 
previous chapter. After non-dimensionalizing all coordinates on the ini- 
tial momentum thickness and the velocity on Au, the mesh size in the cross- 
flow (x^) direction is: ....V 

h 3 = 1.023 

In a mixing layer the eddies are suspected of being elongated in the stream- 
wise direction, so we have set: 


and 

h 2 = h 3 
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We note that if the mixing layer is completely coherent in the spanwise 
direction the size of the mesh in this direction (h^) is not critical. 
The non-dimensional time step was set equal to 


AT = A— 

D. 

in 


0.0799 


6 . 3 Initial Conditions 

We begin by taking the view that the mixing layer is a superposition 
of a random velocity (u,) and a mean velocity profile (u/Au). We want 
the initial random profile to be solenoidal (i.e., V • jj =* 0), random in 
a region of space (see Fig. 6.1), and to decay to zero outside this region. 
In Chapter 3 we showed how to generate an isotropic random velocity 

3 

field ^ on a 16 grid. To generate the random part of the initial 
field that we need here, we start with the field of Chapter 3 and form: 


^(I,J,L) = u x (I,J,L-9) L = 14, ... ,20 


( 6 . 1 ) 


jfe(I,J,L) = 0 


otherwise 


(where I, J, L are the mesh point indices); i.e., a random field over 
the middle of the shear layer that drops abruptly to zero outside. In 
order to smooth out the jump between the two regions, ip is filtered in 
the z-direction with a Gaussian filter. We get: 


i = J i(z’ 


) G(z-z') dz' 


( 6 . 2 ) 


where 


G(z) - 


1 I z 2 \ 
t exp r ~2 ) 

a 3 V 6h~ ' 


The random portion of the initial field is generated by setting 


u = V ip 


(6.3) 
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The initial Conditions were completed by adding to u_ 


an error function 


mean velocity: 


u 

Au 


0.5 erf 


(z-L 3 /2) 
/6 h _ 


Two cases were run: 


Case a: 


max 


Au 


Case b: 


u i 


max 


Au 


= 0.01 (i « 1,2,3) 


= 0.30 (i = 1,2,3) 


(6.4) 


In these two cases the large (grid) structures are assumed to be random 
fluctuations. 

The two-dimensional cases studied in Chapter 5 could be considered as 
unsteady laminar flows, since there is no randomness . We emphasize that 
there are at least two kinds of randomness: 
i) Randomness of the pairing in which the vortices vary in shape, sepa- 
ration distance, strength, number, etc., in a random way. In Chap- 
ter 5 we computed realizations using spacing as the perturbation, 
ii) Randomness meaning noisy (random) fluctuations. 

The calculations described above are designed to look into the second 
type of randomness. To see what the combined effect would be, we ran still 
another case in which the initial field contained a vortex pair (with £ = 
3/16) and a superimposed random field. For the latter, we took the random 
field of case (b) described above. This case will be called (c) . 

Table 6.1 summarizes the cases studied in this chapter. In Appendix D 
we investigate the interaction between streamwise cellular structures and 
spanwise vortex pairing. 


6.4 Momentum Thickness 

In order to study the development of the mixing layer, we would need 
a measure of the effects of the turbulent rotational region on the non- 
turbulent irrotational region; the momentum thickness 0(t) is one such 
measure. We note that Q(t), as defined by Eqn. (5.4) , is a measure of 
the momentum defect of the irrotational region. The momentum defect is due 
to the spreading of vorticity into the irrotational region. Since, in our 
computation, we have dropped the viscous terms, the growth of the momentum 
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thickness measures the inviscid mixing or the entrainment of irrotational 
fluid . 

Figure 6.2 shows the non-dimensional momentum thickness 0/0^ n ^®in 
is the initial momentum thickness) plotted vs. T for the three cases 
considered. We note that in all three cases 0 grows linearly with time. 
The growth rates (d0/Audt) for cases (a) and (b) are not very different, 
despite the large differences in turbulence levels. The values of 0.008 
and 0.011, respectively, are also substantially lower than the growth 
rate (0.018) reported experimentally by S&J; they are, in fact, lower 
than any of the values in Table 1,1. The rate of growth of the momentum 
thickness is only slightly dependent on the intensity of the turbulent 
fluctuations in cases (a) and (b) , and a higher turbulence intensity pro- 
duces a higher growth rate. Furthermore, when large organized structures 
are present (case (c)), the momentum thickness growth rate, d0/Audt = 
0.02 is equal to what it was in the absence of random fluctuations. 

Fig. 6.3 shows the non-dimensional momentum thickness 8/© in plotted vs. 
T, for case (c) and the two-dimeiasional case with 6 = 3/16. Only at the 
early stages of the development of the layer do the random fluctuations 
affect the grovtfth of the momentum thickness . 


6.5 Mean Velocity Profiles 

An important characteristic of the experimental turbulent mixing 
layer is the self-similarity of the mean velocity profiles. In our com- 
putation, the mean velocity < u > is defined by Eqn. (5.8). 

xy 

Figures 6.4a, b, and c show 2 < u > /Au plotted vs. z/0 at 

: xy 

At = 2.4 intervals, for cases (a), (b) , and (c) , respectively. We obtain 
self-similar profiles in all cases. This .means that self- similarity may 
be obtained from a wide variety of different flow structures, and does not 
provide much information about which initial conditions best represent 
physical reality. 


6.6 Mean Turbulence Intensity Profiles 

Experimental observations show that the mean turbulence intensity 
profiles are very nearly self -preserving (Townsend, 1956). This means 
that : 


fc- .. 
£- . 

[• 

? 

i 

*... 

\ ■ 
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= f 


2 (Au) ‘ 


(f) 


(6.5) 


Defining the integral of the turbulent energy 1^ at a given downstream 
distance to be 


-/ 


2 


dz 


- 00 2(Au) 

and substituting (6.5) in (6.6), we get 


( 6 . 6 ) 


h ■ e/ 4 ” '(f) d (f ) - ce 


(6.7) 


where 


C = 


-r 


f(n) d(n) 


Non-dimens ionalizing on the initial integral of the turbulent energy, 


I T,in* weget 


t-t 


T,in 


in 


t. -t 
xn o 


( 6 . 8 ) 


Equation (6.8) shows that I,- grows linearly with time if the profiles of 
2 : 2 ■ 

q /2(Au) are self-similar. To compute 1^, the mean turbulent energy 
defined by Eqn. (5.33) was integrated numerically in the z-direction. 

Figure 6.5 shows I^/l^ plotted vs. T, for the three cases. We 
note that for all three cases I^/l^ decays with time. However, only 
for case (c), in which large structures are present, did the decay level 
off. 

2 2 

Figures 6.6a, b, and c show q /2(Au) plotted vs. z/0, at AT = 

2.4 intervals, for cases (a), (b), and (c) , respectively. Consistent with 
the integral of the turbulence energy results, the turbulence intensity 
decays in time. The most significant drop of the maximum turbulence inten- 
sity occurs in the early stages of the development of the layer. 

The fact that the integral of the turbulence energy decays, instead 
of growing linearly with time, is a clear indication that the term (Eqn. 
(2.16)) used in our equations (2.28) to model the subgrid scale motions. 


55 


.JW 


has too much of an inhibiting effect on the growth of the turbulent fluc- 
tuations. 

In order to support the above argument, we ran a case in which we 
started with the same initial conditions as in case (b), but set C = 0. 

2 2 0 'i 

Fig. 6.7 shows q /2(Au) plotted vs. z/0, and Fig. 6.8 shows I T /I T in 
plotted vs. T, for this case. It is clear that the turbulence intensity 
grows with time, indicating that in case (b) the subgrid scale model is 
inhibiting the growth of the turbulence energy. 

Recall that when the initial conditions contain nothing but large struc- 
tures we obtain self-similar turbulent intensity profiles (see Section 5.8), 
even with C v = 0.188. The decay of the total turbulence energy (Fig. 6.5) 
might suggest that the subgrid scale constant determined for the decay of the 
isotropic turbulence case might be too high for the mixing layer case. How- 
ever, the growth rate of the momentum thickness for case (b) is much lower 
than the growth rate reported experimentally. With = 0, the case (b) 
layer did not grow, i.e., d0/Audt = 0, at least up to T = 9.6, which 
indicates that lowering the subgrid scale constant will not give us a momen- 
tum growth comparable to the experiments. We thus surmise that it is essen- 
tial that large structures be included in the initial conditions if the 
numerical results are to reproduce significant features of the experimental 
mixing layer. In principle, we could begin with a laminar shear layer and 
some small perturbations. The Kelvin-Helmholtz instability would then pro- 
duce large vortical structures and would eventually produce a velocity field 
with the experimentally observed features, A computation of this type would 
require at least an order of magnitude more computing time. As we have 
noted earlier, tue subgrid scale model would inhibit the growth of the per- 
turbations and is not adequate for a computation of transitional flow. We 
shall need to modify the model if transitional flows are to be computed. 

An alternative approach would be to increase the amplitude of the perturba- 
tions and lower the constant of the subgrid scale model, or use a finer mesh. 

6.7 Vorticlty Contours 

In order to investigate the eddy structures and their dynamics, vor- 
ticlty contours in x-z planes have been plotted in Figs. 6.9 and 6.10, for 
the three cases considered, at times T = 0 and T = 16.78. 
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Figure 6.9a shows the spanwise vorticity contours for case (a), at 
time T = 0. The combination of a weak random velocity field and a smooth 
mean velocity distribution yields vorticity contours that are almost unaf- 
; fected by the random fluctuations. The development at T = 16.78, shown 

I in Fig. 6.10a, does not indicate any significant effect of the random 

fluctuations on the mean. The mean field simply masks the weak fluctua- 
i tions in both the initial conditions and at T - 16.78. 

Figures 6.9b show the spanwise vorticity contours for case (b) , at 
I different spanwise (x-z) planes. The combination of a strong random 

velocity field and a mean velocity yields vorticity contours that look 
i spotty. At time T = 16.78, Figs. 6.10b show that the spots appear much 

I 

I more elongated. At some planes (e.g., plane 5), there are two vortex tubes 

! ■ 

! that appear as if they might pair, while other planes show only one vortex 

| tube. This indicates that the initially strong random fluctuations are 

being organized by the mean field, and that the layer is developing through 
a combination of diffusion (due to the subgrid scale model) and vortex pair- 
[.. Tng. 

Figures 6.9c show the spanwise vorticity contours for case (c) at dif- 
ferent spanwise (x-z) planes. Adding random fluctuations to the two span- 
I wise vorticities causes the contour lines of the spanwise vorticity to be- 

come irregular . At time T = 16.78, the vortices have merged in some 
planes (e.g., planes 1-4) in Figs. 6.10a, whereas in other planes (e.g., 
planes 5-6) the vortices are still in the process of merging. This indi- 
cates that strong random fluctuations can affect the dynamics of vortex 
pairing. 



where 


U " = U — < U > 

xy 

Numerically, this quantity is computed as follows. We first calculate 
u", then take its discrete Fourier transform in the y-direction to yield 

VS A . 

u"Cx,k' 2 .»z) • R(x,k 2 ,z) is then defined, to be equal to 

R(x,k2»z) = u"(x,k2,z) u"*(x,k2,z) (6.10) 

A A 

where u" is the complex conjugate of u". Inverse transforming (6.10) 
yields the discrete equivalent of 

R(x,r , z) = f u"(x,y,z) u"(x,y+r,z) dy (6.11) 

Jy 

Finally, line- averaging (6.11) in the x-direction and normalizing yields 
the discrete equivalent to (6.9) . 

Figures 6.11 show R^ at T* = 0 and T = 16.78, plotted vs. r 

at various z locations. We shall define the correlation length to be 

the abscissa of the point where R first crosses the r-axis. 

r uu 

For case (a). Figs. 6.11a show no significant changes in the corre- 
lation length between time T - 0 and T = 16.78. In some parts of the 
flow the correlation length seems to increase, whereas in other parts the 
correlation length seems to decrease. These variations are not signifi- 
cant. 

Figures 6.11b show that when we start with a large random initial 
fluctuation superimposed on a mean profile (case (b)) , the correlation 
length increases with time. This indicates that the layer is becoming 
more organized in the spanwise direction and is consistent with the result 
stated earlier that the vorticity tends to clump. Apparently there is a 
tendency toward the formation of two-dimensional vortices. 

Figures 6.11c show that when we add a random field to coherent struc- 
tures (case (c)), the correlation length decreases slightly with time. The 
only increase in the correlation length occurs at the center of the layer 
(plane 17 in our case) . 
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If the spanwlse correlation length of the streamwise velocity is 
taken as a measure of the coherence of the layer, our results tend to 
indicate that a layer that begins with a random field becomes more co~ 
herent, and one that starts with two-dimensional vortical structures loses 
coherence when the random fluctuations are strong. 

6.9 Summary and Conclusions 

We have shown that the development of the mixing layer is highly de- 
pendent on the initial conditions. This dependence is partly physical and 
partly numerical. Experimentally, the importance of the initial conditions 
on the development of the two-dimensional mixing layer has been pointed 
out by several workers (Bradshaw, 1966; Batt, 1975) . Analytically, the 
subgrid scale models have been developed under the assumption that all the 
energy transferred by the large resolvable scales to the subgrid scales is 
dissipated. The decay of the turbulence intensity in cases (a) and (b) 
indicates that it is doubtful that we can compute transition with the pres- 
ent subgrid scale models. The presence of large structures in the initial 
conditions is essential to the computation of inhomogeneous turbulent flows. 

From the above observations we can conclude that in order to predict 
the initial development of a shear layer one would need a subgrid scale 
model that allows the energy of the small scale field to build up and even- 
tually reach equilibrium with the large eddies. However* the later devel- 
opment of a shear layer can be predicted with the present subgrid scale 
models, provided the large structures are explicitly included in the initial 
conditions. For other flows, it would appear that inclusion of large struc- 
tures that at least approximate those of the physical flow is essential to 
obtain reasonable results. Bass and Orszag (1976) attempted to study the 
evolution of a passive scalar field in a sheared turbulent velocity field, 
but were unable to obtain physically realistic results. This may have been 
due to the omission of the large structures in their initial conditions. 


Chap ter 7 

CONCLUSIONS AND RECOMMENDATIONS 


In this work we have developed an approach to three-dimensional, time- 
dependent computations of flows using the vorticity equations. A general 
method of deriving conservation properties that is applicable to any numer- 
ical method in incompressible fluid mechanics was given; its use simplifies 
the analysis of numerical schemes. 

The use of a filter which is smooth in real space has been shown to 
be essential for the treatment of rotational-irrotational region interac- 
tions. The use o f Fourier tr ansform methods allows accurate and fast treat- 
ment of the term u^UK ” u^j » which arises as a consequence of filtering. 
This is a definite improvement over the expansion in Taylor series (Leonard, 
1973) used in previous studies (Kwak et al., 1975), which we believe should 
be used only when the use of transform methods is not justifiable. 

The vorticity equations have been shown to provide a satisfactory 
basis for the simulation of homogeneous isotropic turbulence. Comparison 
of our results with results obtained using the primitive variable equations 
(Mansour et al. , 1977; Moin et al. , 1978) shows no significant differences. 

A new sub grid scale model has been developed and shown to give results 
comparable to those obtained using the vorticity model (Kwak et al., 1975). 
The new model offers advantages both in computational speed and in storage. 
We found that, for the calculation of isotropic homogeneous turbulence, the 
subgrid scale constant depends only slightly on the numerical method used. 
The variation is about ten percent and is not likely to have a significant 
effect on the computed results in shear flows. The use of Fourier spatial 
differencing has allowed us to look more carefully at the subgrid scale 
model, arid it has been found that replacing exact derivatives with second- 
order differences (roughly equivalent to averaging the model spatially 
(Love and Leslie, 1977)) produces improved behavior of the spectrum. 

No-stress boundary conditions in one direction and periodic boundary 
conditions in the other two directions have been incorporated in a three- 
dimensional, time-dependent code. Flows in which these boundary conditions 
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can be justified (e.g., two-dimensional wakes, planar jets, mixing layer) 
can be investigated using this code. We chose the mixing layer. 

Two-dimensional computations of the turbulent mixing layer have shown 
that pairing vortices produce self-similar mean velocity and turbulence 
intensity profiles. The growth rate of the layer is strongly dependent on 
the initial conditions, a fact also observed experimentally. 

Three-dimensional computations have shown that the presence of large, 
Organized (i.e., not random) structures is essential if the simulation is 
to reproduce the essential features observed in the experiments. These 
computations suggest that in order to simulate the initial development of 
a shear layer one would need a subgrid scale model that allows the energy 
of the small scale field to build up with time and eventually reach equilib- 
rium with the large eddies. However, the later development of a shear layer 
can be predicted with the present subgrid scale models, provided the large 
structures are explicitly included in the initial conditions. 

Our results using different initial conditions indicate that self- 
similarity of the mean velocity profiles can be obtained more easily than 
self-similarity of the turbulence intensities. The addition of strong ran- 
dom fluctuations to a flow containing pairing vortices disturbs the pairing 
in a way that causes the vortex tubes to exhibit spanwise variations, and 
whether or not the merging is completed depends on the spanwise locations. 
This may explain the onset of three-dimensionality seen in experiments. Fig 
7.1 is a conjecture of what we think might happen. The section of the vor- 
tex tube that did not merge could interact with the vortex structure just 
ahead (or just behind) to form a horseshoe vortex. This horseshoe vortex 
may get stretched over several rollers, giving the appearance of cellular 
structures (B&R, Konrad) . 

In Appendix D we study the interaction between streamwise and spanwise 
vorticity. Again, the detailed results depend strongly on the initial con- 
ditions. Each free shear flow is unique, and the universality that is 
sought exists only at large downstream distances. This may mean that the 
computational "prediction'' of free shear flows is feasible only to moderate 
accuracy; the precise behavior of an individual free shear flow may depend 
on physical details that are not easily controlled. This means that some 
experimentation will always be necessary. 
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Work remains to be done on the development of a subgrid scale model 
that incorporates flow-regime dependence. Ideally, one would like a model 
that can handle both transition and developed turbulence. With such a 
model, problems associated with the initial conditions can be studied 
more carefully, since the linear stability theory is well understood and 
the initial conditions can be chosen to be solutions of the Orr-Sotnmerf eld 
equations. This kind of computation will help understand the. effect of the 
initial conditions on the development of the mixing layer, but will not 
reproduce experiments exactly. 

In the case of the mixing layer, the use of periodic boundary condi- 
tions is justifiable only if we move with the mean speed of the flow. How- 
ever, the size of the eddies grows linearly with the streamwise distance 
(in our frame linearly in time), and we reach a point at which the size of 
the box must be increased. In a stationary frame this problem can be 
avoided, but inflow-outflow boundary conditions must be used. We suggest 
that future work should concentrate on developing a method of treating 
the inflow-outflow boundary conditions. 

Eventually, it may be possible to treat practical flows such as air- 
foils, combustion chambers, etc., by these methods. Before that can be 
done, much more effort should first be devoted to developing subgrid scale 
models, treatment of boundary conditions, mesh layout and/or mapping, 
numerical methods, filters, etc., which are the important building blocks 
of large-eddy simulation. 
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Fig. 1.1. r.m.s. streamwise velocity profiles for differ 
ent initial conditions (experimental results 
from Foss, 1977), 
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Fig. 1.2. Mean velocity profiles for different initial 
conditions (experimental results from Foss, 
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Table 4.1 

Computations of the Decay of Isotropic Turbulence 
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Fig. 4.7. Filtered energy spectra. Pseudo-spectral corapu- 
■ ■ tation with 32 mesh; 2nd-otder differencing 

' model w-2. C = 0.188. 
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Fig. 5.3. Computational box and coordinate system. 
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Fig. 5.4. Initial conditions setup |B 
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Fig. 5.5. Non-dimensional momentum thickness (0/0 in) as a function 
of time for various 3. Two-dimensional computations. 




Fig. 5.6a. Contour plots of the spanwise vorticity (aig) 
Cor 3 = 6/16, at time T = 0. Constant vor- 
ticity lines ate plotted at eight levels. 
Higher numbers on these Lines indicate higher 

levels (w_ = 0.702). 
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Contour plots of the spanwise vorticity (o^) 
3 = 4/16, at time T = 0. Constant vorticit 
lines are plotted at eight levels. Higher ni 
on these lines indicate higher levels (oi„ 
0.444). 2 > ms 
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Fig. 5.6d, Contour plots of the. s panwise vorticity (ta*) for 

3 = 3/16, at time T — 0. Constant vorticity lines 

are plotted at eight levels. Higher numbers on these 

lines indicate higher vorticity levels (oi rt - 0.421). 
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Fig. 5. 6e. Contour plots of the spanwise vorticity ((tig) for 

3 = 2/16, at time T - 0. Constant vorticity lines 

are plotted at eight levels. Higher riumbers on 

these lines indicate higher levels (w 0 „ = 0.416) 

° Z.max 


















Fig. 5.7a. Contour plots of the spanw ise vorticity (toj) for 

3 = 6/16, at time T = 16.78. Constant vorticity 

lines are plotted at eight levels. Higher numbers 

on these lines indicate higher levels (to,, 
n on/, \ Z, max 
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Fig. 5.7b. Contour plots of the spanwise vorticity (o^) for 
3 = 5/16, at time T = 16.78. Constant vorticity 
lines are plotted at eight levels. Higher numbers 
on these lines indicate higher levels (w ? - 

0.358). 




Fig. 5. 7c 


Contour plots of the sp.an.wise vortic i ty (w?) Tor 
3 = 4/16, at time T - 16.78. Cons taut vorticity 
lines are plotted at eight levels. Higher numbers 
on these lines indicate higher vorticity levels 
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Fig. 5.7d. Contour plots of the spanwise vorticity (0)2) for 
g = 3/16, at time T = 16.78. Constant vorticity 
lines are plotted at eight levels. Higher numbers 
on these lines indicate higher levels = 

0.276). 




i 

i; 



Fig. 5 . 7 o . Contour plots of the spanwlse vorticity ( 0 J 2 ) 

for (5 = 2/16, at time T = 16.78. Constant 

vorticity lines are plotted at eight levels. 

Higher numbers on these lines indicate higher. 

levels (to. = 0,248) . 
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Fig. 5.7f . Contour plots of the spnnwise vortieity (t0 9 ) for 
P » 1/16, at; ttm' ?•* 16.78, Constant vorticity 
lines are plotted at eight: levels. Higher numbers 
On these lines indicate higher levels (Wi 'm 
0.245) . 2,raax 
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Fig. 5.8. Mean velocity profiles. Ttoo-d imens ional computation 
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Fig. 6.1. Three-dimensional computation box. Random velocity 
setup and coordinate system. 
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Table 6.1 


Three-Dimensional Computations of Turbulent Mixing Layers 


Case 

Amplitude of 

Random Field 

Initial Conditions 

a 


u f 
X 

max 

= 0.01 

Random field + mean 

Au 

b 


u i 

max 

» 0.30 

Random field -f mean 

Au 



u i 

max 

» 0.30 

Random field + 2 spanwise 
vortices (3 = 3/16) 

c 

Au 



Fig. 6.2. Non-dimensional momentum thickness (0/0^) ns a func- 
tion of time. Three-dimensional computations. 
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Fig. 6.4a. Mean velocity profiles, 
(case a) . 
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Fig. 6.4b. Mean velocity profj 
(case b) . 
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Fig. 6.6b. Mean turbulence intensi' 
computations (case b). 
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Fig. 6.7. Mean turbulence intensity profiles. Three-dimensional 
computation (C^ = 0). 
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Pig. 6.8. Integral of the turbulence energy as a function of 
time. Three-dimensional computation (C =0). 




Fig. 6.9a. 


Contour plots of the spanwise vorticity (u^) in 
an x-z plane, at time T = 0. Constant vorticity 
lines are plotted at eight levels. Higher numbers 
on these lines indicate higher vorticity levels 
^2, max = °- 228 > case a ) • 
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Figs. 6.9b. Contour plots of the spanwise vorticity ( 0 ) 3 ) for different x-z 


planes, at time T = 0. In each plane, constant vorticity lines 
are plotted at eight levels. Higher numbers on these lines in- 
dicate higher vorticity levels (case b) . 
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14.322 ; al 2 - 0.427 = 15.345 
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Figs. 6.9b (continued) 














Figs. 6.9c. Contour plots of the spanwise vorticity (a) 2 ) for different x-z 
planes, at time T = 0. In each plane, constant vorticity lines 
are plotted at eight levels. Higher numbers on these lines indi- 
cate higher vorticity levels (case c) . 
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Figs. 6.9c (continued) 
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Figs. 6.9c (continued) 










Fig. 6.10a. Contour plots of the s panwise vorticity (tog) in an 
x-z plane, at time T = 16.78. Constant vorticity 
lines are plotted at eight levels. Higher numbers on 
these lines indicate higher vorticity levels 

(W 2,max = °- 185 ’ case a) ‘ 
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Figs. 6.10b (continued) 
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Figs. 6.10b (continued) 
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mixing layer, Az/6. ^ = 1.023; case 
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Fig. 7.1. Formation of streamwise cellular structures in 
a mixing layer . 
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Appendix A 


SUBGRID SCALE MODELS FOR THE VORTICITY EQUATIONS 

In Chapter 2 we propose to use the following models for W^. (Eqn. 
(2.9)): 


Model 0)-l 


Model w-2 


where 


w u - 

“ £ ijk ax^ (2V T S k^ 

(2.15) 

W . . = - 

3x . (V T W i^ + 3x. (V T W j ) 
j i 

(2.16) 

V T 

= (C/) 2 (to^)* 5 

(2.14) 

ii 

1 / a - . a - \ 

2 \8x. u i 8x. u j / 

(2.12) 


The models should satisfy the following necessary conditions: 

1. they should be antisymmetric, 

2. they should vanish in an irrotational region, and 

3. they should be an energy sink. 

Condition 1 is readily seen to be satisfied by these models. We note 
that in an irrotational region, ok = 0. Hence, = (C^A) (mho^) 2 - 0, 
and the model vanishes in an irrotational region; i.e., condition 2 is also 
satisfied . 

In order to show the dissipative nature of the subgrid scale models 
W-l and to-2, consider the following equation: 


3io. 
x 

at 


a 


ax. 

i 



(A.l) 


where the nonlinear terms in Eqn. (2.28) have been dropped. Multiplying 
Eqn. (A.l) by ip , and integrating over the flow volume, we get: 
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f* iH“i dv ■ -fh 

We want to show that Eqn. (A. 2) reduces to 


. ~ — W. . dv 

1 dX. XJ 

j 


3 n — 
t y 2 i a 


J 

at 


dv 


- e 


(A. 2) 


(A. 3) 


where e 0 . 

Model t)-l 

Substituting Eqn. (2.15) in Eqn. (A. 2) for W _ , integrating by parts, 
and using periodic boundary conditions, we obtain 

i dv ‘ “ 2 / v i dv (A - 4) 

and we have for this case: 

e - 2 / v dv i 0 

since >_ 0 , 


Model W-2 


In a similar way, substituting Eqn. (2.16) into Eqn. (A. 2) for W^_. , 
we can show that 

u.u. dv = - / \> 0).0). dv (A. 5) 

dt^/ 2 x x ^ T x x 

and we have for this case 




w.O). dv > 0 

xx — 
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Appendix B 


Fast Discrete Sine Transform (FDST) 

The discrete analogs to the expansion in Fourier sine series (Eqns. 
(3.14) and (3.15)) are 


f(x) - 


f (n) sin 


• r n7rx ~ l 

ln L(N-l)hJ 


f (n) = 


(N-l) 4. 

3=o 


S f(x) sln 


where n 
h 
x 
N 
L 


o,i,. v , N-i, 

L/ (N-l) , 

jh, j = 0,1,... , N-l, 

number of mesh points, 

length of the computational box. 


(3.16) 


(3.17) 


Both the forward and inverse sine transforms involve identical sums. 
Eqn. (3.] 7) can be rewritten as: 

~2(N-1)-1 


“ CN^lT Im 


^ F(x) exp: 

L J =0 


( ~2t r inx \ ' 
2(N-l)h J 


(B.l) 


where 


F(x) = f(x) j = 0,1,..., N-l, 

= 0 j = N, N+l, . . . , 2(N-1)-1 


We note that the summation 

2 (N -l 1 -1 

4 

3 


2 F(x) exp ( icS-ijh) 


(B.2) 


is equivalent to (3.12) with N^ = 2(N-1) , and an FFT routine can be 
used to evaluate this sum. 
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Fast Discrete Cosine Transform (FDCT) 


The discrete analog to the expansion in Fourier cosine series (Eqns. 
(3.21) and (3.22)) are: 


f(x) = ^ f?(n) cos («J?L5h) 


n=o 


(3.19) 


where 


aC 

f’(n) 


f’(x) = 


f (n) = 


fef (n) 
t f C (n) 
(hi (x) 
t f(x) 


( S £,(x) cos [wm) (3 - 20) 


(N-l) 


n = 0,N-1 
n = 1, . . , ,N-2 
j = 0,N-1 
j = l,...,N-2 


where 


n = 0, ...,N-1 

h ^ L/(N-1), 

x - jh j = 0, . . . ,N-1, 

N = number of mesh points, 

L — length of the computational box. 

Both the forward and inverse transforms involve identical sums, 
Eqn. (3.19) can be rewritten as: 

~2(E-!)-l 


f (x) - - Re 


11=0 


„ / s /— 2 tt inx \ 

,W exp ( 2(N-l)h / 


(B • 3) 


where 


F(n) = H f (n) n = 0,N-1, 

-.V'- ' f C (n) n = 1. . . . ,N-2, 

■ 0 ; n = N, . . . ,2(N-1)-1. . 

We note that the sum in (B. 3) is identical to the sum(B. 2), and an 
FFT routine can be used to evaluate it. In fact, the sine and cosine 
transforms can be done simultaneously, if it is necessary to have both. 
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Appendix C 


Effect of a Sinusoidal Vorticity Perturbation on a Uniform Vortex Array 

In Chapter 5 we have studied the effect of perturbing a uniform array 
of vortices by offsetting the spacing of the vortices (3 > 0). In this 
appendix we study the effect of adding a sinusoidal vorticity perturbation 
to a uniform array of vortices ( 3 = 0 ). 

1. Initial Conditions 

The initial conditions studied in this appendix were generated by 
starting with a uniform array of point vortices on the centerline of our 
computational box: 

oo 2u = C^ ^SCx-L^/4) + SCx-SL^M)^ 6(z~L 3 /2) (C.l) 

We then add a cosine vorticity distribution to (C.l): 

0 ) £ = W 2u ” C 2 cos ^ C z— L 3 / 2 ) (C.2) 

Eqn. (C.2) is then filtered with a relatively wide Gaussian filter (Eqn. 
(5.19)) to yield the initial conditions. The initial velocity is then 
non-dimens ionalized on Au and the length scales on 0^. The computa- 
tional details, i. e. , number of mesh points, mesh size, time steps, and 
boundary conditions, are the same as in Chapter 5. Only the initial con- 
ditions were changed. 

2. Results . 

The momentum thickness ( 6 ) is defined by Eqn. (5.4). Fig. C.l shows 
0/6 in plotted vs. T for ^ C 1 = 0.1/20, 1/20, 2/20, 4/20. We note that 
the growth rate of the layer is highly dependent bn the strength of the 
perturbation. The growth rate more than doubles from 0.016 to 0.035 when 
the strength of the perturbation is doubled (^/Cj from 2/20 to 4/20). 
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We note also that for high amplitude perturbations, C 2^ C ± = the 

growth rate starts to level off for T > 12.0. This saturation is also 
observed experimentally by Oster et al. (1978); they have oscillated the 
initial conditions of a two-dimensional mixing layer. 

Figures C.2 and C.3 show the non-dimensional mean velocity and turbu- 
lence intensity (as in Sections 5.7 and 5.8) plotted vs. z/0 for C 2^ C 1 = 
2/20. We note that the mean velocity profiles are self-similar. This is 
not surprising, since self-similarity of the mean velocity profiles is 
easily obtained (see Section 6.5). Turbulence intensity profiles (Fig. C.3) 
show that self-similarity is also more or less obtained for the present 
case. 

These results are similar to those obtained in Chapter 5 by using a 
spacing perturbation. Apparently the perturbation can take any of a number 
of forms, and the characteristics of the shear layer will be nearly the 
same. Under experimental conditions, the nature of the perturbation is 
difficult to determine. What we do note is that reproduction of the ex- 
perimentally observed growth rate does require large perturbations, which 
are apparently created by either the inflow or outflow conditions of the 
experiment. 



Appendix D 


INTERACTIONS BETWEEN STREAMWISE AND SPANWISE VORTICITY 

In Chapter 6 we studied the effect of a random fluctuation on vortex 
pairing. In this appendix we study the interactions between a streamwise 
cellular vortex structure and spanwise vortex pairing. 


1. Initial Conditions 

The initial conditions studied in this appendix were generated by 
adding to a row of spanwise vortices (3 = 3/16) a row of streamwise vor- 
tices of alternating signs: 


= C 2 sin 


expl- 


(z-iy 2) 


(D.l) 


The same computational setup described in Chapter 6 is used, i.e., the same 
boundary conditions, number of mesh points, mesh sizes, and time step. 

Figure D.l shows a contour map in the y-z plane of the streamwise 
vorticity. We note that 00^ displays a cellular structure and that 
does not initially have a streamwise variation. We ran two cases: 


Case a: 


Case b ; 


1 max 

2 max 


1 max 

2 max 


= 0.037 


= 0.370 


2. Results 

We first look at the development of the momentum thickness, 0(t), de- 
fined by Eqn. (5.4), in time. The non-dimensional mean velocity (Section 
5.7) and mean turbulence intensity (Section 5.8) are also considered. The 
interaction between the spanwise vortices and the streamwise vortices is 
studied using contour plots. Note that we have a three-dimensional box 
and that contour plots in different planes for different vorticity direc- 
tions will be considered. 


Figure D.2 shows 0/0^ n plotted vs. T. The momentum thickness 
growth rate, d0/&udt = 0.020, for Case (a) is the same as it was in the 
absence of the streamwise vortices. However, the momentum thickness 
growth rate, d0/Audt - 0.040, doubled for Case (b) , 

Figures D.3a and -b show 2< u > /Au plotted vs. z/0 for Cases (a) 
and (b) , respectively, at At = 2.4 intervals. We note that both cases 

produce self-similar mean velocity profiles. 

2 2 

Figures D.4a and -b show q /2(Au) plotted vs. z/0 for Cases (a) 
and (b) , respectively, at AT = 2.4 intervals. The mean turbulence inten- 
sity results for Case (a) are similar to those we obtained when the stream- 
wise vortices were not present. As in the 2-D case (with 0 = 3/16), the 
mean turbulence intensity decays slightly, then reaches a self-similar 
situation. For Case (b) , in which we have strong streamwise vortices, 

Fig. D.4b shows that the turbulence intensity grows with time, and the pro- 
files do not show self-similarity. 

(a) Contour Plots in the x-z Planes 

Figures D.5 show constant vorticity contours of the spanwise ( 0 ) 2 ) 
vorticity at time T = 16.78. In both cases the spanwise vortices have 
paired. The shapes are similar, but the roller is slightly distorted for 
Case (b) as compared to Case (a) and to the 2-D results (see Fig. 5,7d). 
This indicates that the streamwise vortices did not affect the merging of 
the spanwise vortices, but the strong streamwise vortices (Case (b)) have 
affected the shape of the roller. 

Figures D.6 show constant vorticity contours of the streamwise vortic- 
ity for Cases (a) and (b) . These figures indicate that the streamwise vor- 
tices have been convected to the edges of the mixing layer by the spanwise 
vortices. There is also clear evidence of vortex stretching. 

Figure D.7 shows the projection of the vorticity vector at T = 16. 78-, 
for Case (b) . We can see clearly that the originally straight vortex lines 
have been convected and stretched by the spanwise vortices to assume an 
inverted S shape. 
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Figure D.8 shows constant vorticity contours of the spanwise vortic- 
ity for Case (b) . The spanwise vortices have been convected and stretched 
by the strong counter-rotating streamwise vortices and exhibit spanwise 
waviness. This means that the contact area between the rotational fluid 
and the irrotational fluid has increased, which leads to an increase in 
the entrainment rate. This waviness also explains the increase in the 
turbulence intensity and high growth of the momentum thickness of the mix- 
ing layer. J te that the mean quantities are defined as horizontal planar 
averages and, with this definition, the wavy layer appears thicker and more 
turbulent than a strictly two-dimensional layer. 

The above results indicate that the effect of the streamwise vorticity 
on the spanwise vorticity is almost independent of the effect of the span- 
wise vorticity on the streamwise vorticity. Indeed, a straight line of 
particles placed at the center of the layer in the streamwise direction 
would be convected to form an inverted S shape in the presence of the two- 
dimensional vortex pairing. A straight line of particles initially passing 
through the center of an array of counter-rotating vortices will be con- 
vected to assume a wavy shape. 


i 



1.6 


1.5 


1.4 



0 2.4 4.8 7.2 9.6 2 12.0 14.4 16.8 


. T 

Fig. C.l. Non-dimensional momentum thickness (0/0. ) as a func- 
tion of time for various C^/C^. 
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Fig. D.5a. Contour plots of the spanwise vorticity (w 2 ) in an x-z 
plane (y/0 in = 4.09), at time T = 16.78. Constant vor 
ticity lines are plotted at eight levels. Higher numbers 
on these lines indicate higher levels (case a). 



Fig. D.5b. Contour plots of the spanwise vorticity (c^) in an x-z 
plane (y/®i n ~ 4.09), at time T = 16.78. Constant vor 
ticity lines are plotted at eight levels. Higher numbers 
on these lines indicate higher levels (case b) . 




Fig. D.6a. Contour plots of the streamwise vorticity (to-t) in an x 
plane (y/®i n = 4*09), at time T — 16.78. Constant vor 
ticity lines are plotted at eight levels. Higher numbers 
on these lines indicate higher vorticity levels (case a) . 



Fig. D.6b. Contour plots of the streamwise vorticity (w^) in an 
x-z plane (y/9 in = 4.09), at time T = 16.78. Con- 
stant vorticity lines are plotted at eight levels. 
Higher numbers on these lines indicate higher levels 
(case b). 
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Fig. D.7. Projection of the vorticity vector in an x-z plane 
(y/^in = 4*09) at time T = 16.78. 
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Appendix E 


c COMPUTER PROGRAM WRITTEN TO CALCULATE TURBULENT MIXING LAYERS 

C 

C 

KCOMDECK AVG 

COMMON/AVG/ AVG1 , AVG2 , AVG3 , CCF 
KCOMDECK BLANK 

COMMON DUDX(16, 16,33) 

KCOMDECK DATA7 

C0MM0N/DATA7/ FR< 16 , 16 ) , FI ( 16 , 16 ) 

KCOMDECK DATA9 

C0MM0N/DATA9/ IMAX , JMAX, LMAX 
KCOMDECK DAT21 

C0MM0N/DAT21/ XR(6A) ,XI(64) 

KCOMDECK DEL 

COMMON/DEL/ DELTAX, DELTAY, DELTAZ 
KCOMDECK DIM 

COMMON/DIM/N1 ,N2 , N3 
KCOMDECK FLT 

COMMON/FLT/ FILT1 <16 ) , FI LT2 ( 16 ) , FILT3< 33 ) 

KCOMDECK LARGE2 

COMMON/ LARGE2/ U( 16 , 16 , 33 ) , V ( 16 , 16 , 33 ) , W( 16 , 16 , 33) 

LEVEL 2,U, V,W 
KCOMDECK LARGE3 

COMMON/L ARGE3/ GU( 16 , 16 , 33 ) , GV ( 16 , 16 , 33 ) , GW( 16 , 16 , 33 ) 

LEVEL 2 , GU , GV , GW 
KCOMDECK LARGES 

COMMON/L A RGE5/ 01 ( 16 , 16 , 33 ) , 02C 16 , 16 , 33) , 03< 16 , 16 , 33 ) 

LEVEL 2,01,02,03 
KCOMDECK MEANVOR 

COMMON/MEANVOR/ V0R(32,33) 
xrnMnFrx pp 

COMMON/PR/ CCPW,CCPF,CCPD 
KCOMDECK WV 

COMMON/WV/ WAVEXC 16 ) , WAVEYC 16 ) > WAVEZ( 33 ) , WAVEXS( 16 ) , WAVEYS ( 16 ) 

1 , WAVEZSC 33 ) 

KCOMDECK XL 

COMMON/XL/ XPART (160),YPART(160),ZPART(160), NCHAR( 160) 

KDECK MAIN 

PROGRAM MA I N(INPUT, OUTPUT, TAP E8,TAPE9,TAPE10) 
CXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXKXXKXKXKXKXKXKXX 


C MAIN CONTROLS THE COMPUTATION SEQUENCES. K 
C IN THIS ROUTINE WE ADVANCE IN TIME . K 
C THE EXTERNALS USED IN THIS ROUTINE ARE K 
C CFILTER K 
C CONVEC K 
C CURLO K 
C DATARED K 
C EDVIS K 
C INVERS K 
C MEANINI K 
C MOVLEV . ■ X 
C SFILTER K 
C SFILTER K 
C SGS K 
C STFILT K 
C STREAD K 
C STWV K 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

INTEGER TIME, TSTART, TEND 
COMMON/TIM/ TSTART, TEND 

COMMON/LARGE4/ RU ( 1 6 , 16 , 33 ) , RV< 16 , 16 , 33 ) , RWC 16 , 16 , 33) 

COMMON/NORM/ DELU, THETA * 

LEVEL 2 , RU , RV , RW 
COMMON/ DA TCNT/ IDATCNT 
KCALL MEANVOR 
KCALL XL 

COMMON/CONST/C1 0 0,C101,IJK,IJ ,NHP1 , HALF 
KCALL DAT21 
KCALL LARGE2 
KCALL BLANK 


REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 


XCALL DIM 
XCALL AVG 

C START THE READOUT OF INPUT 
CALL STREAD 

C SET THE COEFICIENT OF THE SUBGRID SCALE MODEL 

C — g 

C SEt'cOF = 1 FOR THE FIRST TIME STEP 
COF=1.0 
IJ=N1XN2 
IJK=N1XN2XN3 
DO 1 L=1,LMAX 
DO 1 J = 1 , JMAX 
DO 1 1 = 1 , IMAX 
U(I,J,L>=0. 

VCI,J,L)=0. 

W(I,J,L>=0. 

01 ( I » J > L ) = 0 . 

02(1, J,L)=0 . 

03( I > J , L ) =0 . 

1 CONTINUE 

C SET THE WAVE NUMBERS 
CALL STWV 

CXXXXXSET THE INITIAL CONDITIONS 
CALL MEANINI 

C NON DIMENSIONALIZE THE TIME STEP ON DELU/THETA 
DT=0.03125*DELU/ THETA 
C SET THE FILTER WIDTH = 2XMESH SIZE 
AVG1=2 . 0 
AVG2=2 . 0 
AVG3=2.0 

C0EF2=(CX(AVGlXDELTAXxAVG2XDELTAYXAVG3XDELTAZ)xx(l./3. )>XX2 
DO 123 L=1,LMAX 
DO 123 J=1 , JMAX 
DO 123 1 = 1, IMAX 
RU( I , J , L ) =0 . 

RV(I,J,L)=0. 

RW( I , J , L ) = 0 . 

123 CONTINUE 
ICOUNT=0 
TIME=0 

C WRITE ON TAPE 9 TO BE STORED ON DISC PACK 
PRINT 1100 ,TIME 

WRITEC9) TIME,01,02,03,DT,DELTAX,DELTAY,DELTAZ,DELU,THETA 
IDATCNT=0 

C COMPUTE THE STATISTICS OF THE INITIAL CONDITIONS 
CALL DATARED 

C SET THE FOURIER TRANSFORM OF THE GAUSSIAN FILTER 
CALL STFILT 
IDUM=30 

DO 300 TIME =T ST ART, TEND 

CXXXXXCOMPUTE THE ADVECTIVE AND STRETCHING TERMS 
CALL CONVEC 

CALL SFILTER(GU, DUDX, N1 , N2 , N3) 

CALL S FI LT ER( GV , DUDX , N1 , N2 , N3 ) 

CALL CFI LTERC GW , DUDX , N1 , N2 , N3 ) 

CXXXXXCOMPUTE THE EDDY VISCOSITY 

CALL EDVIS ( C0EF2 , DUDX , N1 , N2 , N3 ) 

CXXXXXCOMPUTE THE SGS MODEL 

CALL SGS ( U , V , W, N1 , N2 , N3 ) 

CXXXXXADVANCE IN TIME 
DO 80 0 L = 1 , LMAX 
DO 800 J=l, JMAX 
DO 800 1=1, IMAX 


XCALL LARGE5 
XCALL DEL 
XCALL DATA9 
XCALL WV 
XCALL DATA7 
XCALL FLT 
XCALL LARGE3 
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01CI,J,L)=01(I, J,L)+DTX(COFxGU(I,J,L)“0.5XRU(I, J,L)) 
02CI,J,L)=02CI, J,LMDTX(CQFXGVCI,J,L)-0.5XRVCI,J,L)) 
03(I,J,L)=03CI,J,L)+DTX(COFXGWCI,J,L)-0.5XRW(I,J,D) 

800 CONTINUE 

CXXXXXSTORE THE PREVIOUS TIME STEP 

CALL M0VLEV(GU(1,1»1),RU(1,1»1),IJK) 

CALL MOVLEV ( GV ( 1 , 1 , 1 ) , RV ( 1 , 1 ,1 ) , I JK) 

CALL MOVL EV(GW( 1,1,1), RW< 1,1,1), I JK) 

CXXXXXTHE VORTICITY AT THE NEXT TIME STEP HAS BEEN COMPUTED 
CxxxxxFIND THE CORESPONDING VELOCITY FIELD 
CALL INVERS(01,GU,DUDX,1,N1,N2,N3) 

CALL INVERSC02,GV,DUDX,2,N1,N2,N3) 

CALL IN VERS (03,GW,DUDX,3,N1,N2,N3) 

CALL CURL0<GU,GV,GW,U,V,W,N1,N2,N3) 

C SET COF = 1.5 FOR SUBSEQUENT TIMES ( ADAMS-BASHFORTH) 

C0F=1 . 5 

ICOUNT=ICOUNT+1 
IICOUNT=ICOUNT-IDUM 
IF CIICOUNT .NE. 0) GO TO 300 
ICOUNT=0 
PRINT 1100, TIME 
WRITE( 9 ) TIME, 01, 02, 03 
CALL DATARED 
300 CONTINUE 

1100 FORMAT ( 1H1 , 5X, X TIME STEP =X, 15) 

1000 FORMAT (1P8E15.7) 

STOP 

END 

XDECK CFILTER 

SUBROUTINE CFILTER (HR, HI, N1,N2,N3) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C CFILTER COMPUTES THE FILTER OF THE HR VARIABLE BY EXPANDING IN 

C FOURIER SERIES IN THE X- AND ¥- DIRECTIONS AND FOURIER COSINE 

C SERIES IN THE Z-DIRECTION 

C THIS ROUTINE USES AS EXTERNALS 

C FDOT 

C FFTX 

C FFTY 

C A CALL TO STFILT INITIATE THE VALUES OF FILT1 , FILT2, AND FILT3 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXX 

DIMENSION HR(N1,N2»N3),HI(N1,N2»N3) 

XCALL DATA9 
XCALL FLT 
XCALL DATA7 
XCALL DAT21 

LEVEL 2, HR 

CC=1 . 0/( IMAXXJMAX) 

IJ=N1XN2 
DO 10 J = 1 , UMAX 
DO 10 1=1, IMAX 
DO 20 L=1 , LMAX 
XRCL)=HR(I, J,L) 

20 CONTINUE 

CALL FDCT (1.0) 

DO 30 L = 1 , LMAX 
HI ( I , J , L ) =XR< L ) 

30 CONTINUE 
10 CONTINUE 

DO 90 L = 1 , LMAX 

CALL MOVLEV(HI(l»I, L ) , FR ( 1 , 1 ) , I J ) 

CALL FFTX(l.O) 

CALL FFTY(1. 0,1.0) 

DO 50 J=1 , JMAX 
DO 50 1=1, IMAX 

FRC I, J ) =FR ( I > J ) XFltTl ( I ) XFILT2( J ) XFILT3( L ) 

FI C I » J ) =FI( I , J JXFILTI ( I )XFILT2( J )XFILT3( L ) 

50 CONTINUE 

CALL FFTXI-l .0) 

CALL FFTY(-1 . 0 , CC) 

CALL MOVLEV(FR(l,l),HI(l,l,L),IJ) 


40 CONTINUE 

DO 60 J=1,JMAX 
DO 60 1 = 1 , IMAX 
DO 70 L = 1 , LMAX 
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DO 60 1=1, IMAX 
DO 70 1=1, IMAX 
XR(L)=HI(I»J»L) 

70 CONTINUE 

CALL FDCT(-l.O) 

DO 80 1=1, IMAX 
HR<I, J,L)=XR(L) 

80 CONTINUE 
60 CONTINUE 
RETURN 
END 

XDECK CONVEC 

SUBROUTINE CONVEC 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXX 
C THIS SUBROUTINE COMPUTES THE CONVECTIVE AND STRETCHING 
C TERMS AND STORES THEM IN GU,GV,GW 

C THIS ROUTINE USES AS EXTERNALS 

C COSPART 

C PARTIAL 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

XCALL LARGE2 

XCALL LARGE3 

XCALL LARGE5 

XCALL BLANK 

XCALL DATA9 

XCALL DIM 

IJK=N1XN2XN3 

CXXXXXTERM FOR THE X-DIRECTION 
DO 10 L=1,LMAX 
DO 10 J=1,JMAX 
DO 10 1=1, IMAX 

GUM, J,L)=U<I, J,L)X02(I,J,L)-VCI,J,L)X0lU, J,L) 
GVCI,J,L)=U<I,J,L)X03(I,J,L)-W(I, J,L)X01(I, J,L) 

10 CONTINUE 

CALL PARTI AL ( 2 , GU, Nl , N2, N3 ) 

CALL M0VLEV(DUDX(1,1,1),GU(1,1,1),IJK) 

CALL C0SPART(GV,N1,N2,N3) 

DO 20 L = 1 , LMAX 
DO 20 J=1 , JMAX 
DO 20 1=1, IMAX 

GUC I , J , L ) =GU( I , J , L )+DUDX( I , J , L ) 

20 CONTINUE 

CXXXXXTERM FOR THE Y-DIRECTION 
DO 30 L = 1 , LMAX 
DO 30 J = 1 , JMAX 
DO 30 1=1, IMAX 

GV(I, J,L)=V(I, J,L)X01(I, J,L)-UCI, J,L)X02(I, J,L) 

GW(I, J,L)=V(I, J,L)X03(I, J,L)-W(I, J,L)X02(I, J,L) 

30 CONTINUE 

CALL PARTIAL<1,GV,N1,N2,N3) 

CALL MOVLEV(DUDX(l,l,l),GV( 1,1,1), IJK) 

CALL COSPART C GW, N1,N2,N3) 

DO AO L=l, LMAX 
DO 40 J=l, JMAX 
DO 40 1=1, IMAX 

GV ( I , J , L ) =GV ( I , J , L ) +DUDXC I , J , L ) 

40 CONTINUE 

CXXXXXTERM FOR THE Z-DIRECTION 
DO 50 L=1 , LMAX 

■ DO 50 J = l, JMAX . ■ ■ . . 

DO 50 1=1, IMAX 

GW (I , J , L )=W( I , J , L )X01 C I , J , L )-U(I , J , L )X03C I , J , L) 

UCI, J,L)=W(I, J,L)X02(I, J,L)-V(I,J,L)X03(I, J,L) 

50 CONTINUE 

CALL PARTIAL! 1 , GW, N1 , N2 , N3 ) 
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CALL MOVLEV(DUDX(l,l,l),GW(ld,l),IJK) 

CALL PARTIAL(2,U,N1 ,N2,N3) 

DO 60 L=1,LMAX 
DO 60 J=1.JMAX 
DO 60 1=1 , IMAX 

GWd , J > L ) =GWC I , J , L ) + DUDX( I , J , L ) 

60 CONTINUE 
RETURN 
END 

XDECK COSPART 

SUBROUTINE COSPART( U , N1 , N2 , N3 > 

CXXXXXXXXXXKXXXXXXXXXXXXXXXXXXHXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C COSPART COMPUTES THE PARTIAL IN THE Z-DIRECTION OF U BY EXPANDING 
C IN FOURIER COSINE SERIES 

C THE EXTENALS USED IN THIS SUBROUTINE ARE 

C FDCT 

C FDST 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DIMENSION U(N1,N2,N3) 

XCALL BLANK 
XCALL DAT21 
XCALL WV 
XCALL DATA9 

LEVEL 2 » U 
DO 10 J=1 , JMAX 
DO 10 1=1, IMAX 
DO 20 L=1 , LMAX 
XR( L ) =U ( I , J , L ) 

20 CONTINUE 
SIGN=1 . 0 
CALL FDCT ( SIGN ) 

DO 30 L=1 , LMAX 
XR(L)=-XR(L)XUAVEZ(L) 

30 CONTINUE 
SIGN=-1.0 
CALL FDST (SIGN ) 

DO AO L = 1 » LMAX 
DUDX(I, J,L)=XR(L) 

AO CONTINUE 
10 CONTINUE 
RETURN 
END 

XDECK CURLO 

SUBROUTINE CURLO ( U , V , W, 01 , 02 , 03 , N1 , N2 , N3 ) 

DIMENSION 01 ( N1 , N2 , N3 ) , 02 ( N1 » N2 , N3 ) , 03 (Nlj N2 , N3 > 

DIMENSION U(N1,N2,N3),V(N1,N2,N3),W(N1,N2,N3) 

XCALL DATA9 
XCALL BLANK 

LEVEL 2,U,V,W,01 ,02,03 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXX 
C THIS SUBROUTINE COMPUTES THE CURL OF THE VORTICITY FIELD 

C THE EXTERNALS USED IN THIS ROUTINE ARE 

C PARTIAL 

C SINPART 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
CXXXXXCURL IN THE X-DIRECTION 
IJK=N1XN2XN3 

CALL PART I AL ( 2 , W, N1 » N2 , N3 ) 

CALL MO VL EV( DU DX (1,1,1), 01(1, 1,1),IJK) 

CALL SINPART(V,N1,N2,N3) 

DO 10 L=1 , LMAX 
DO 10 J=1 , JMAX 
DO 10 1=1, IMAX 

01(1, J,L)=0 1(1, J,L)~ DUDX(I, J,L) 

10 CONTINUE 

CXXXXXCURL IN YHE Y-DIRECTION 
CALL SINPART(U,N1,N2,N3) 

CALL MOVLEV( DUDX( 1,1,1), 02(1, 1,1), IJK) 

CALL PARTI AL (1 ,W, N1 , N2 ,N3) 

DO 20 L=1 , LMAX 




DO 20 J=1,JMAX 
DO 20 1 = 1 , IMAX 

02(1, J,L) =02(1, J,L)-DUDX< I, J,L) 

20 CONTINUE 

CXXXXXCURL IN THE Z-DIRECTION 

CALI PARTIAL(1»V»N1»N2,N3) 

CALL MOVLEV ( DUDX( 1,1,1), 03( 1,1,1). I JK) 

CALL PARTIAL(2,U,N1,N2,N3> 

DO 30 L = 1 » LMAX 
DO 30 J=1,JMAX 
DO 30 1=1, IMAX 

03(I,J,L)=03(I»J»L)-DUDX(I,J|L) 

30 CONTINUE 
RETURN 
END 

XDECK CURLU 

SUBROUTINE CURLU( U , V , W, 01 , 02 , 03 , N1 , N2 , N3 ) 

DIMENSION 01(N1,N2,N3> , 02( N1 , N2 , N3 ) , 03 ( N1 , N2 , N3 ) 

DIMENSION U(N1,N2,N3>,V(N1,N2,N3),W(NI,N2,N3) 

XCALL DATA9 
XCALL BLANK 

LEVEL 2,U, V,W,01,02,03 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


C THIS ROUTINE COMPUTES THE CURL OF THE VELOCITY FIELD X 
C AND STORES IT IN 01,02,03, X 
C THE EXYERNALS USED IN THIS SUBROUTINE ARE X 
C COSPART X 
C PARTIAL x 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

CXXXXXCURL IN THE X-DIRECTION 
IJK=N1XN2XN3 

CALL PARTIAL(2,W,N1,N2,N3) 

CALL MOVL EV( DUDX( 1 » 1 , 1 ) > 01 (1* 1 , 1 ) > I JK ) 

CALL COSPART (V, N1 , N2 , N3 ) 

DO 10 L = 1 , LMAX 
DO 10 J=1 , JMAX 
DO 10 1=1 , IMAX 

01(I,J,L)=01CI,J,L>- DUDX(I, J,L) 

10 CONTINUE 

CXXXXXCURL IN YHE Y-DIRECTION 
CALL C0SPART(U,N1,N2,N3) 

CALL MOVL EV ( DUDX( 1, 1,1), 02(1,1,1), I JK) 

CALL PARTI AL < 1 , W, N1 , N2 , N3 ) 

DO 20 L=1 , LMAX 
DO 20 J = 1 , JMAX 
DO 20 1=1, IMAX 

02(1, J,L)=02(I, J , D-DUDXCI, J,L) 

20 CONTINUE 

CXXXXXCURL IN THE Z-DIRECTION 

CALL PARTIAL ( 1 , V,N1 , N2 , N3 ) 

CALL MOVLEV(DUDX( 1,1,1), 03(1, 1,1), IJK) 

CALL PARTIAL(2,U,N1,N2,N3) 

DO 30 L=1 , LMAX 
DO 30 J=1 , JMAX 
DO 30 1=1, IMAX 

03(1, J,L)=03(I,J,L)-DUDX(I, J,L) 

30 CONTINUE 
RETURN 
END 

XDECK DATARED 

SUBROUTINE DATARED 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C THIS SUBROUTINE COMPUTES THE STATISTICS OF THE COMPUTATION 


C USUM = PLANAR AVERAGE OF THE STREAMWISE VELOCITY X 
C VSUM=PLANAR AVERAGE OF THE SPANWISE VELOCITY X 
C WSUM = PLANAR AVERAGE OF THE CROSSFLOW VELOCITY x 
C 01 SUM = PLANAR AVERAGE OF THE STREAMWISE VORTICITY X 
C 02SUM = PLANAR AVERAGE OF THE SPANWISE VORTICITY X 
C 03SUM = PLANAR AVERAGE OF THE CROSSFLOW VORTICITY X 
C USQ = R.M.S STREAMWISE VELOCITY X 
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C VSQ = R.M.5. OF THE SPANWISE VELOCITY ETC... 

C UVSTRE5 = PLANAR AVERAGE OF U3W3 

C PLOVALE = VOLUME AVERAGE OF THE TOTAL ENERGY 

C ENERGY = INTEGRAL OF THE TURBULENCE ENERGY 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

XCALL DEL 
XCALL MEANVOR 
XCALL PR 
XCALL LARGE2 
XCALL LARGE3 
XCALL LARGE5 
XCALL DAT21 
XCALL DATA9 
XCALL BLANK 
XCALL DIM 

DIMENSION USM(33) , V5M(33) ,U5(33) »V5(33) » WS ( 33 )>01S( 33) >025(33) 

1 , 035(33), ES(33),ENS<33),Z0(33) 

DIMENSION DUMSPC33) 

COMMON/D ATCNT/ IDATCNT 

IDATCNT=IDATCNT+1 

LMAXM1=LMAX-1 

C3=l ,/LMAXMI 

CN0RM3=1./(IMAXXJMAX) 

PRINT 1100 
UTOT=0. 

VTOT=0. 

WTOT=0. 

OITOT-O. 

O2TOT=0. 

O3TOT=0. 

OVRALE=0. 

TOTENER=0. 

TOTENST-O . 

DO 100 L=1,LMAX 
USUM=0. 

VSUM=0. 

WSUM=0. 

01 SUMF O'. 

023UM=0. 

03SUM=C. 

USQ = 0 . 

VSQ=0. 

WSQ=0. 

015(5 = 0. 

02S<5 = 0. 

03SQ=0 . 

ENERGY=0 . 

ENSTRCP=0 . 

UVSTRES=0 . 

PLOVALE=0 . 

DO 110 J = l, JMAX 
DO 110 1=1, IMAX 
USUM=U5UM4-U(I, J,L) 

VSUM=VSUM+V(I, J,L) 

WSUM=WSUM+W( I , J , L ) 

01SUM=01SUMf01(I, J,L) 

02SUM=02SUM+02(I,J,L) 

03SUM=03SUM+03 ( I , J , L ) 

110 CONTINUE 

USUM=U5UMXCN0RM3 
VSUM=V5UMXCN0RM3 
WS U M = W 5 UM x C N 0 RM3 
01SUM=0ISUM*CN0RM3 
02SUM=02SUMXCNORM3 
03SUM=035UMXCN0RM3 
GW( 1 , 1 , L ) =USUM 
GW(2, 1 , L )=VSUM 
GWC 3 , 1 , L ) =WSUM 
DO 160 J=l, JMAX 
DO 160 1=1, IMAX 
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xx*::*; 


USQ=USQ+CU(I, J,L)-USUM>K*2 
VSQ=VSQ+(V(I, J,L)-VSUM)X*2 
WSQ=WSQ+ (W( I * J > L )~WSUM ) HX2 
01SQ=01SQ+<01<I, J,L)-01SUM)X*2 
025Q = 02SQ4-(Q2(I,J>L ) “02 SUM)* #2 
03SQ = 03SQ+(03(I»J»L)“03SUM)XX2 
UVSTRES=UVSTRES+(UCI,J,L)-USUM>X(W(I, J,L)-WSUM> 
P10VALE=P10VALE+(U(I, J,L)**2+V(I, J,L)X*2+U(I, J,L)*X2) 
160 CONTINUE 

USQ=USQXCN0RM3 

VSQ=VSQ*CN0RM3 

WSQ=WSQ*CN0RM3 

01SQ=01SQXCN0RM3 

02SQ=02SQ*CN0RM3 

03SQ=03SQ*CN0RM3 

ENERGY = ( USQ+VSQ+WSQ )X0 . 5 

ENSTRQP=C01SQ+02SQ+03SQ)*0.5 

UVSTRES=UVSTRES*CN0RM3 

USQ=SQRT ( USQ ) 

VSQ=SQRTCVSQ) 

W5Q=5QRT tWSQ ) 

01SQ=SQRT C 01 SQ ) 

02SQ=SQRT(02SQ) 

03SQ=SQRT(03SQ) 

US ( L ) =USQ 
VSCL)=VSQ 
WS(L)=WSQ 
015 (L)=01SQ 
Q2S ( L ) =02SQ 
03S(L)=03SQ 
USMCL)=USUM 
VSM<L)=VSUM 
ESC l ) =ENERGY 
ENS(L)=ENSTROP 
GWC A » 1 , L ) =01 SUM 
GWC 5 > 1 » L ) =02SUM 
GWC6, 1,L)=03SUM 
GVC1 t 1 t L) = USQ 
GV C 2 , 1 , L ) =WSQ 
XI C L ) =UVSTRES 
CC“1. 

IFCl’.EQ. 1) CC = 0 . 5 
IFCL .EQ. LMAX) CC=0.5 
OVRALE=OVRALE+PLOVALE*CCX0.5 
UTOT-UTOT+USUMXCC 
VT OT = VTOT + VSUflxcC 
WTOT-WTOT+WSUM*CC 
01T0T=01T0T +01SUMXCC 
02T0T=02T0T+02SUM*CC 
03T0T=03T0T+03SUM*CC 
TOTENER-TOTENER+ENERGYXCC 
TOTENST=TOTENST+ENSTROPXCC 
100 CONTINUE 

UT0T=UT0T#C3 

VT0T=VT0T*C3 

WT0T=WT0T*C3 

01T0T=01T0T*C3 

02T0T =02T0T*C3 

03T0T =03T0TXC3 

TQTENER=T0TENER*C3 

T0TENST=T0TENST*C3 

DELU=GWC 1,1* LMAX) -GWC 1,1,1) 

DELU=1./DELU 

THETA=(0.25-(GWC1,1,1)XDELU)XX2)X0.5 
DO 170 L = 2 , LMAXM1 

THETA=THETA+(0.25-(GW(l,l,L)XDELU)XX2) 

170 CONTINUE 

THETA=THETA+CO .25-CGWC1 ,1,LMAX)XDELU)XX2)X0 .5 

THETA=THETAXDELTAZ 

DO 300 L r l , LMAX 
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ZOC L ) =( L-C ( LMAX-1 )>2U ) ) XDELTAZ/THETA 
USMCL )=USMCL)XDELUX2.0 
VSMCL)=VSMCL)XDELU 

GW( 5, 1 , L ) =GW< 5 » 1 , L )*DELU*THETA 

XI ( L ) "XI C L ) x ( DELUXX2 ) 

USCL)=USCL)XDELU 

VS(L)=VSCL)XDELU 

WSCL)=WSCL)XDEIU 

03SU>=03SCL)x(DELU*THETA> 

01S(L)=01S(L)X(DELUkTHETA) 

02SC L ) =02SC L )XC DELUXTHETA ) 

ESC L ) =ES C L ) XDELUXX2 
ENSCL)=ENSCL)XCDELUXTHETA)XX2 

PRINT 3000,USMCL),VSMCL),XICL),GWC5,1,L),USCL),VS(L)>WS(L), 
1 015CL),02SCL),03SCL),E5(L),ENS(L),Z0(L) 

300 CONTINUE 

WRITEC8 ) USM,VSM,XI,U5,VS,WS,01S,02S,03S,ES, ENS»ZO» THETA 
PRINT 1700, THETA 
PRINT 1200 

PRINT 1000,UTOT,VTOT,WTOT,OlTOT,O2TOT,O3TOT,TOTENER,TQTENST 
PRINT 2A00 ,OVRALE 

2400 FORMAT C IX, X OVER ALL ENERGY IN COMPUTATION BOX =X,1PE15.7> 
DO 180 L = 1 , LMAX 
DO 180 1=1 , IMAX 
GU(I,1,L)=0. 

GU(I,2,L)=0. 

GUCI,3,L)=0. 

GUCI,4,L)=0. 

180 CONTINUE 

C10Y=1 ./FLOAT C JMAX) 

DO 190 L = 1 > LMAX 
DO 190 J = 1 , JMAX 
DO 190 1=1, IMAX 

GUC I , 1 , L ) =GU (I,1,L)+02(I,J,L)*C10Y 
GUCI,2,L)=GUCI,2,L)+UCI, J,L)XC10Y 
GUC 1 , 3 » L ) =GU C 1 , 3 , L )+W( I , J , L ) XC10Y 
190 CONTINUE 

DO 230 L = 1 , LMAX 
DO 230 1=1, IMAX 
GUC I » 2 > L ) =GUC 1 , 2 , L )-GWC 1 , 1 , L ) 

GUC 1 , 3 , L ) =GUC I,3,L)-GW(3,1,L) 

230 CONTINUE 
PRINT 2200 

2200 FORMAT C 1H1 , IX, X LINE AVERAGE OF VORTICITY*) 

PRINT 2300,CCCGUCI,1,L),I= 1,16 ) , L ) , L =1 , LMAX) 

IFCCCPD . N E . 1 . ) GO TO 240 
PRINT 2500 

2500 FORMAT C 1H1 , IX, X LINE AVERAGE OF U-COMPONENT X)- 
PRINT 2300 , C C CGUC 1 ,2, L) , 1 = 1,16 ) , L ) , L = 1 , LMAX ) 

PRINT 2500 

PRINT 23 0 0 , CC CGUC 1,2, L) , I = 17, IMAX) ,L),L = 1, LMAX) 

PRINT 2600 

2600 FORMAT C 1H1 , IX, X LINE AVERAGE OF W-COMPONENT x) 

PRINT 2300,(CCGUCI,3,L),I= 1,16 ) , L ) , L-l , LMAX) 

PRINT 2600 

PRINT 2300, CCCGUCI, 3, L), 1=17, IMAX), L),L=1, LMAX) 

2300 FORMAT C1X,16F8.3,I3) 

240 CONTINUE 
PRINT 2000 
DO 250 L =14,20 
DO 260 1=1, IMAX 
XRC I > -GUC I 
XICI)=0. 

260 CONTINUE 

CALL FFTCXR,XI,IMAX,-1) 

IFC IDATCNT . EQ . 1 ) DUM5PC L)=SQRTCXRC2)XX2+XIC2)XX2) 
IFCDUMSPCL) .LT. 0.0000001) GO TO 250 
DO 270 1=1, IMAX 

XRCI)=S<3RTCXRCI)XX2+XICI)XX2)/DUMSPCL) 

270 CONTINUE 
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PRINT 1800,L»(XR(I),I = 1>8) 

250 CONTINUE 

1800 FORMAT ( 1X>XWV02X, 15, 1P8E14.6) 

1C00 FORMAT (IP8E15.7) 

1100 FORMAT (2X,*USUM*,6X,*VSUMk,5X,XUWSTR#,5X,X02SUMX,7X,XUSQ#,7X,*VSQ# 

1,7X,XWSQX,6X,X01SQX,6X»X02SQX,6X»X03SQX > 5X,XENERGYX,AX, XENSTROPX , 3 
2X,*PIANE*) 

1200 F0RMAT(///,1X,X UTOT IN X-Y VTOT IN X-Y WTOT IN X-Y Ol TOT 

1 IN X-Y 02TOT IN X-Y 03TOT IN X-Y X) 

1300 FORMAK1X,# USUN IN Y-Z VSUM IN Y-Z WSUM IN Y-Z OlSUM IN Y 

1-Z 02SUM IN Y-Z 035UM IN Y-Z X) 

1A00 FORMAT!/// , IX , X UTOT IN Y-Z VTOT IN Y-Z WTOT IN Y-Z OITOT 

1 IN Y-Z 02T0T IN Y-Z 03T0T IN Y-Z KJ. 

1500 FORMAT ( IX, x USUN IN Z-X VSUM IN Z-X WSUM IN Z-X OISUM IN Z 

1-X 02SUM IN Z-X 03SUM IN Z-X X) 

1600 FORMAT (/// , IX, X UTOT IN Z-X VTOT IN Z-X WTOT IN Z-X OITOT 

1 IN Z-X 02T0T IN Z-X 03TOT IN Z-X X) 

1700 FORMAT ( IX, X MOMENTUM THICKNESS X,1PE15.7> 

2000 FORMAT ( 1H1 ) • 

3000 FORMAT (1P13E10.2) 

C TEST THE DIVERGENCE OF THE VELOCITY AND VORTICITY FIELDS 
CALL DIV 

C TEST THE SOLUTION OF THE POISSON EQUATIONS ,I.E. THAT X 

C THE CURL OF OUR VELOCITY FIELD IS EQUAL TO THE VORTICITY FIELD 

CALL CURLU<U,V,W,GU,GV,GW,N1,N2,N3) 

ERRMAX1 =0 . 

ERRMAX2=0. 

ERRMAX3=0 . 

DO 17 L = 1 , LMAX 
DO 17 J = 1 , UMAX 
DO 17 1 = 1 , IMAX 

GU(I,J,L)=ABS(GU(I,J,L)-01(I,J,L)) 

GV<I, J,L)=ABS(GV<I, J,L)-02CI, J,L)> 

GW( I » J » L ) =ABS(GW( I , J , L )-03( I , J » L ) ) 

IF (GUCI,J,L) .GT. ERRMAX1 ) ERRMAX1 =GU( I , J , L ) 

IF CGVCI,J,L) .GT. ERRMAX2) ERRMAX2=GV ( I , J , L ) 

IF ( GW( I > J , L ) .GT. ERRMAX3 ) ERRMAX3=GW( I , J , L ) 

17 CONTINUE 

PRINT 1110 , ERRMAX1 , ERRMAX2, ERRMAX3 

1110 FORMAT ( IX, X ERRMAX1 »x , E15 . 7 , XERRMAX2 =X,E15.7,X ERRMAX3 =X,E15.7) 
RETURN 
END 

XDECK DIV 

SUBROUTINE DIV 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C THIS ROUTINE TESTS THE DIVERGENCE OF THE VELOCITY FIELD M 

C AND THE DIVERGENCE OF THE VORTICITY FIELD . X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

XCALL LARGE2 

XCALL LARGE3 

XCALL LARGE5 

XCALL BLANK 

XCALL DATA9 

XCALL DIM 

IJK=N1XN2XN3 

CALL PARTIAL ( 1 , U ,N1 , N2,N3 ) 

CALL MOVL EV( DUDX( 1 , 1 , 1 ) , GU( 1 , 1 , 1 ) , I JK) 

CALL PARTIAL (2 , V, N1 , N2 , N3 ) 

CALL MOVLEVC DUDX( 1 , 1 , 1 ) , GV ( 1 , 1 , 1 ) , I JK ) 

CALL S INPART (W, Nl, N2,N3) 

DIVMAX=0. 

DO 1 L = 1 , LMAX 
DO 1 J=1,JMAX 
DO 1 1=1, IMAX 

DUM-ABS C GU ( I, J , L ) +GV( I , J , L >+DUDX( I , J , L ) ) 

IF (DUM .GT. DIVMAX) DIVMAX=DUM 
1 CONTINUE 

PRINT 1100, DIVMAX 

CALL PARTIAL ( 1 ,01 , Nl , N2 , N3) 

CALL MOVL EV ( DUDX( 1 , 1 , 1 > , GU ( 1 » 1 , 1 ) , I JK) 
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CALL PARTIAL(2,02,N1 »N2 » N3 ) 

CAL L MOVL EV ( DUDX<1 , 1 , 1 ) , GV C 1 , 1 , 1 ) , I JK ) 
CALL C0SPART(03,N1,N2,N3) 

DIVMAX-0. 


DO 2 L = 1 > LMAX 
DO 2 J=1,JMAX 
DO 2 1 = 1 » IMAX 

DUM=ABS(GU (I»J,L)+GV(I,J,L ) + DUDX( I , J , L ) ) 

IF (DUM .GT. DIVflAX) DIVMAX=DUM 
2 CONTINUE 

PRINT 1200 , DIVMAX 
1100 FORMAT ( IX ,x MAXIMUM VELOCITY 
1200 FORMAT ( IX, X MAXIMUM VORTICITY 
RETURN 
END 

KDECK EDVIS 

SUBROUTINE EDVIS ( COEF2 »E»N1 »N2 ,N3) 

DIMENSION E(N1,N2,N3) 

c X X X X X X X X XXXX X X XXX X X X X XX X X XXXXX X XXXX XXXX X XXXX XXXX X XX XX X X X X X X X X X X X XX X X X XX 

: THIS SUBROUTINE COMPUTES THE EDDY VISCOSITY AND STORES IT IN E * 


DIVERGENCE =*,E15.7> 
DIVERGENCE =*,E15.7> 


CXX XXXX XXX XXXXXX XXX XXX XXX XXX XXX XXX XXXX XXX XXX XXXXXXXXX XXX XXX XXXXXXXXXXXXX 

*CALL LARGE5 
XCALL DATA9 

DO 3 L = 1 » LMAX 
DO 3 J=1,JMAX 
DO 3 1=1, IMAX 

E(I,J>L) = 01(I,J,L)XX2+02(I»J,L)XX2+03(I»J»L)XX2 
E(I, J,L)=SQRT(ECI,J,L)>*C0EF2 
3 CONTINUE 
! RETURN 
END 

#DECK FDCT 

SUBROUTINE FDCT (SIGN) 


CKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XXX XXXXXXXXXXXXXX XXXXX XX 


C FDCT COMPUTES THE FAST DISCRETE COSINE TRANSFORM OF THE VARIABLE 
C XR AND STORES IN XR 


CXXXXX XXXX XXXXXXX XXXXX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XXX 


XCALL DAT21 
XCALL DATA9 

LM1 =LMAX~1 
CC=2./FL0ATCLM1) 

LL=2*LM1 

XR(l)=XR(l)/2. 

XR( LMAX) =XR( LMAXJ/2 . 

LP1=LMAX+1 
DO 1 L=LP1,LL 
XR(L>=0. 

1 CONTINUE 
DO 3 L=1,LL 
XI ( L ) =0 . 

3 CONTINUE 
ISN--SIGN 

CALL FFT(XR , XI , LL , ISN ) 

IF (SIGN .GT. 0.) GO TO 200 
DO 100 L = 1 » LMAX 
XR(L)=XR(L)XCC 
100 CONTINUE 
200 CONTINUE 
RETURN 
END 

xnprK Fn«iT 

SUBROUTINE FDST(SIGN) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


C FDST COMPUTES THE FAST DISCRETE SINE TRANSFORM OF THE VARIABLE * 
C XR AND STORES IT IN XR . . , * 

CXX XXXXX XXXXXXX XXX XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XX XXXXX XXXXX XXX 


XCALL DAT21 
XCALL DATA9 

LM1=LMAX-1 

CC=2./FL0AT(LM1) 
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LL=2XLM1 
XRC1>=0. 

XR(LMAX>=0 . 

LP1=LMAX+1 
DO I L = LP1 > LL 
XR< L ) =0 . 

1 CONTINUE 

DO 3 L = 1 , L L 
XI ( L ) =0 . 

3 CONTINUE 
ISN=-SIGN 

CALL FFT(XR,XI,LL,ISN) 

IF (SIGN .GT. 0. ) GO TO 200 
DO 100 L =1 , LMAX 
XKL >=XI(L)XCC 
100 CONTINUE 
200 CONTINUE 

DO 2 L =1 » LMAX 
XR(L)=-SIGNXXI(L) 

2 CONTINUE 
RETURN 
END 

XDECK FFT 

IDENT FFT <A,B,N,ISN) 

ENTRY FFT 

x RADIX 2 COMPLEX FAST FOURIER TRANSFORM, COMPUTED IN PLACE. 

* SEE ON COMPUTING THE FAST FOURIER TRANSFORM, R. SINGLETON, 
x COMM. ACM, V.10, N.IO, PP.6A7-65A, OCT. 1967. 

X ARRAY A CONTAINS THE REAL COMPONENT OF THE DATA AND RESULT, 

X ARRAY B CONTAINS THE IMAGINARY COMPONENT. 

X N, THE NUMBER OF COMPLEX DATA VALUES, 

X MUST BE A POWER OF 2 AND GREATER THAN 1 

x THE SIGN OF ISN IS THE SIGN OF THE EXPONENTIAL IN THE TRANSFORM. 
X THE MAGNITUDE OF ISN IS THE INCREMENT SIZE FOR INDEXING 
X A AND B, AND IS ONE IN THE USUAL CASE. 

X DATA MAY ALTERNATIVELY BE STORED FORTRAN COMPLEX 
X IN A SINGLE ARRAY, IN WHICH CASE THE MAGNITUDE 

x OF ISN IS TWO AND ADDRESS B IS A(2), I.E. 

X CALL FFT2( A , A ( 2 ) , N , 2 ) 

x INSTEAD OF 


FFT2C 2 
FFT2C 3 
FFT2C A 
FFT2C 5 
FFT2C 6 
FFT2C 7 
FFT2C 8 
FFT2C 9 
FFT2C 10 
FFT2C 11 
FFT2C 
FFT2C 
FFT2C 
FFT2C 
FFT2C 
FFT2C 
FFT2C 


X 


CALL 

FFT2C A , B , N , 1 ) 


FFT2C 

X 

PROGRAM CONTAINS SINE TABLE FOR 

MAXIMUM N OF 32768 

FFT2C 

X 

6A00 

TIME 

FOR 

N-1.02A, 220 M.SEC. 

FFT2C 

X 

6 A00 

TIME 

FOR 

N=2**M IS 21.5XNXM MICRO-SEC. 

FFT2C 

X 

6600 

TIME 

FOR 

N=102A, AA M.SEC. 


FFT2C 

X 

6600 

TIME 

FOR 

N=2XXM IS A.3XNXM 

MICRO-SEC. 

FFT2C 

X 

RMS 

ERROR 

FOR 

TRANSFORM-INVERSE 

IS LESS THAN 1.3E-13 

FFT2C 

X 

FOR N=32768 

OR SMALLER. 


FFT2C 

X 

FORTRAN 2. 

3 SUBROUTINE 


FFT2C 

X 

BY R 

. C. SINGLETON, STANFORD RESEARCH INSTITUTE, NOV. 1968 

FFT2C 

L100 

SXO 


B3 

NN 

FFT2C 



SB A 


BO 

KK = 0 

FFT2C 



SB3 


B3-B7 

NN=NN-INC 

FFT2C 



AXO 


1 

KSPAN=NN/2 

FFT2C 



SB5 


BO 

K2 = 0 

FFT2C 



SB6 


XO 


FFT2C 



SX1 


B5 

K2=K2 

FFT2C 



EQ 


B6 , B7 , FFT 

I F< KSPAN .EQ. INC) RETURN 

FFT2C 

L110 

SBA 


B3-BA 

KK=NN-KX 

FFT2C 



SB5 


B3-B5 

K2=NN-K2 

FFT2C 



SA2 


Bl + BA 

EXCHANGE A(KK),A(K2) AND 3(KK) ,B(K2) 

FFT2C 



SA3 


B1 + B5 


FFT2C 



SAA 


B2 + BA 


FFT2C 



NX7 


X2 


FFT2C 



SA5 


B 2 + B 5 


FFT2C 



NX6 


X3 


FFT2C 



SA7 


A3 


FFT2C 



SA6 


A2 


FFT2C 



NX7 


XA 


FFT2C 



NX6 


X5 


FFT2C 


12 
13 
1 A 

15 

16 

17 

18 
19 


21 
22 
23 
2 A 

25 

26 

27 

28 

29 

30 

31 

32 

33 
3A 

35 

36 

37 

38 

39 
AO 
A 1 
A2 
A3 
AA 
A5 
A6 
A7 
AS 
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L120 


1130 


FFT 


L10 


L20 


130 




i*v' Y,j ff ,a 


'JSSL 


SA7 

A5 



FFT2C 

A9 

SA6 

AA 

END OF EXCHANGE 


FFT2C 

50 

IT 

B6 , BA , L 1 10 

IFCKSPAN .LT. KK) GO 

TO L110 

FFT2C 

51 

SBA 

BA+B7 

KK=KK+INC 


FFT2C 

52 

SB5 

B6 + B5 

K2=KSPAN+K2 


FFT2C 

53 

SA2 

B1 + B4 

EXCHANGE ACKK),A<K2) 

AND B(KK) »B(K2) 

FFT2C 

5A 

SA3 

B1 + B5 



FFT2C 

55 

SAA 

B2 + BA 



FFT2C 

56 

NX7 

X2 



FFT2C 

57 

SA5 

B2+B5 



FFT2C 

58 

NX6 

X3 



FFT2C 

59 

SA7 

A3 



FFT2C 

60 

SA6 

■A2 



FFT2C 

61 

NX7 

XA 



FFT2C 

62 

SXO 

B6 

K=KSPAN 


FFT2C 

63 

HX6 

X5 



FFT2C 

6 A 

SA7 

A5 



FFT2C 

65 

SA6 

AA 

END OF EXCHANGE 


FFT2C 

66 

AXO 

1 

K = K/2 


FFT2C 

67 

1X1 

Xl-XO 

K2=K2-K 


FFT2C 

68 

PL 

XI , 1130 

IF(K2 .GE. 0) GO TO 

L130 

FFT2C 

69 

LXO 

1 

K = K+K 


FFT2C 

70 

SB* 

B4 + B7 

KK=KK+INC 


FFT2C 

71 

1X1 

Xl+XO 

K2=K2+K 


FFT2C 

72 

SB5 

XI 

K2=K2 


FFT2C 

73 

GE 

B5 t8 A r L110 

IF(K2 .GE. KK) GO TO 

L110 

FFT2C 

7 A 

LT 

B4,B6,L120 

IF< KK .LT. KSPAN) GO 

TO L120 

FFT2C 

75 





FFT2C 

76 

SB1 

XI 



INSR1 

1 

SA1 

Al + 1 



INSR1 

2 

5B2 

XI 



INSR1 

3 

SA1 

Al + 1 



IN5R1 

A 

SB3 

XI 



INSR1 

5 

SA1 

Al + 1 



INSR1 

6 

SB4 

XI 



INSR1 

7 


SAA 

BA 

ISN 

FFT2C 77 

MX2 

1 

MASK 

FFT2C 78 

SA5 

L60 


FFT2C 79 

SA3 

B3 

N 

FFT2C 80 

LX2 

57 


FFT2C 81 

PX7 

X3 


FFT2C 82 

BX6 

-X2XX5 


FFT2C 83 

PL 

XA,L10 

IFCISN .GE. 0) GO TO L10 

FFT2C 8 A 

BX6 

X2+X5 


FFT2C 85 

BXA 

-XA 

INC--INC 

FFT2C 86 

LX3 

32 


FFT2C 87 

SA6 

A5 


FFT2C 88 

NXO 

B5,X3 


FFT2C 89 

PX2 

XA 


FFT2C 90 

SB7 

XA 


FFT2C 91 

DX7 

X2XX7 


FFT2C 92 

SA1 

B5+S 

S(M) 

FFT2C 93 

SB3 

X7 

NN=INCXN 

FFT2C 9A 

SB6 

X7 

KSPAN=NN 

FFT2C 95 

EQ 

LAO 

GO TO LAO 

FFT2C 96 

SA3 

CD 


FFT2C 97 

RXA 

X2*X1 

SDXCN 

FFT2C 98 

RX7 

X2XX0 

SDXSN 

FFT2C 99 

RX5 

X3*X0 

CD*SN 

FFT2C100 

RX6 

X3XX1 

CDXCN 

FFT2C1 0 1 

RXA 

XA-X5 


FFT2C102 

RX6 

X6+X7 


FFT2C1 03 

NX5 

XA 


FFT2C1 0 A 

RX7 

X1-X6 


FFT2C105 

RXO 

X0+X5 


FFT2C106 

NXl 

X7 


FFT2C107 

SB5 

B6 + BA 

K2=KSPAN+KK 

FFT2C108 

SA2 

Bl + BA 

A(KK) 

- FFT2C109 

SA3 

B1 + B5 

A C K2 ) 

FFT2C110 

SAA 

B2 + BA 

B(KK) 

FFT2C111 
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RX6 

X2+X3 

I 

FFT2C1 12 

SA5 

B2 + B5 

B(K2) 

FFT2C113 

RX2 

X2-X3 

RE . : 

FFT2C11A 

SA6 

A2 

ACKK? 

i 

t 

FFT2C115 

RX7 

X4+X5 

FFT2C116 

RX3 

X1*X2 

CNXRE 

FFT2C1 1 7 

RX4 

XA-X5 

IM 

FFT2C118 

SA7 

AA 

B(KK > 

FFT2C1 1 9 

RX5 

XOXXA 

SN*IM 

FFT2C120 

RX2 

X0XX2 

SNXRE 

FFT2C121 

RX6 

X3-X5 


FFT2C122 

RXA 

X1*XA 

CN*IM 

FFT2C1 23 

SA6 

A3 

ACK2 ) 

FFT2C12A 

RX7 

X2+X4 


FFT2C125 

SBA 

B6 + B5 

KK=KSPAN+K2 

FFT2C126 

SA7 

A5 

B(K2) 

FFT2C127 

LT 

BA > B3> L30 

IF<KK .LT. NN) GO TO L30 

FFT2C128 

SB5 

BA-B3 

K2=KK-NN 

FFT2C129 

BX1 

-XI 

CN=-CN 

FFT2C130 

SB A 

B6-B5 

KK=KSPAN-K2 

FFT2C131 

LT 

B5,B4,L30 

IFCK2 .LT. KK) GO TO L30 

FFT2C1 32 

SBA 

BA+B7 

KK=KK+INC 

FFT2C133 

SA2 

SD 


FFT2C134 

LT 

B4,B5,L20 

IFCKK .IT. K2) GO TO L20 

FFT2C1 35 

SBA 

BO 

KK=0 

FFT2C136 

SX5 

B6 


FFT2C137 

AX5 

1 

KSPAN=KSPAN/2 

FFT2C138 

SB6 

X5 


FFT2C139 

SB5 

B6 + BA 

K2=KSPAN+KK 

FFT2C1A0 

SA2 

Bl + BA 

ACKK) 

FFT2C1A1 

SA3 

B1 + B5 

ACK2) 

FFT2C1A2 

SAA 

B2 + BA 

BCKK) 

FFT2C1A3 

RX6 

X2+X3 


FFT2C1AA 

SA5 

B2 + B5 

B(K2) 

FFT2C1A5 

RX7 

X2-X3 


FFT2C1A6 

SA6 

A2 

ACKK) 

FFT2C1A7 

SA7 

A3 

ACK2) 

FFT2C1A8 

RX6 

XA+X5 


FFT2C149 

SBA 

B6 + B5 

KK=KSPAN+K2 

FFT2C150 

RX7 

XA-X5 


FFT2C151 

SA6 

AA 

BCKK) 

FFT2C152 

SA7 

A5 

BCK2) 

FFT2C153 

LT 

B A > B 3 > L 5 0 

IFCKK .LT. NN) GO TO L50 

FFT2C154 

EQ 

B6 >B7, L100 

IFCKSPAN .EQ. INC) GO TO 

L100 FFT2C1 55 

SA1 

A1 

SCM) 

FFT2C156 

SBA 

B7 

KK=INC 

FFT2C157 

RX6 

X1XX1 


FFT2C158 

SA1 

Al + 1 

M=M+1, SCM) 

FFT2C159 

FX6 

X6+X6 


FFT2C160 

SA3 

ONE 


FFT2C16 1 

SA6 

CD 

CD-2*S CM ) XX2 

FFT2C162 

BXO 

XI - v 

SN=SD 

FFT2C163 

RX6 

X3-X6 

CN=1.0-CD 

FFT2C16 A 

BX7 

xo 


FFT2C165 

NX1 

X6 


FFT2C166 

SA7 

SD 


FFT2C167 

EQ 

L30 

GO TO L30 

FFT2C168 

DATA 

9. 587379909597 7 3 A6E-5 

FFT2C169 

DATA 

1 . 917 A7597 31 070331 E-A 

FFT2C17 0 

DATA 

3.83A9518757 1 39559E-A 

FFT2C171 

DATA 

7 . 66 9903187 A270A53E-A 

FFT2C172 

DATA 

1 . 533980 18628 A76 56 E~3 

FFT2C1 73 

DATA 

3. 067956762 96 59763 E- 3 

FFT2C17A 

DATA 

6.13588A6A915AA75AE-3 

FFT2C175 

DATA 

1 .2271538285719926E-2 

FFT2C176 

DATA 

2 . A5A1228522912288E-2 

FFT2C177 

DATA 

A.906767A327A1801AE-2 

FFT2C178 

DATA 

9.80171A0329560602E-2 

FFT2C179 

DATA 

1.95090 32201612827 E-l 

FFT2C180 

DATA 

3.82683A3236508977E-1 

FFT2C181 
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DATA 

0.7071067811865475 

FFT2C182 

ONE 

DATA 

1.0 

FFT2C183 

CD 



FFT2C184 

SD 

END 


F FT 2 C 1 8 5 

XDECK 

FFTX 

SUBROUTINE 

FFTX(SIGN) 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C FAST FOURIER TRANSFORM IN X-DIRECTION 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XCALl DAT A9 
XCALL DAT A7 
XCALL DAT21 

ISN=-SIGN 

IF (SIGN .LT. 0.) GO TO 3 
DO 1 J=l, JMAX 
DO 1 1 = 1 , IMAX 
FI(I, J)=0. 

1 CONTINUE 
3 CONTINUE 

DO 100 J=1 , JMAX 
DO 110 1=1, IMAX 
XR( I )=FR( I , J ) 

XI ( I )=FI ( I , J ) 

110 CONTINUE 

CALL FFTCXR, XI , IMAX, ISN ) 

DO 120 1=1, IMAX 
FR( I , J )=XR( I ) 

FICI, J)=XI(I) 

120 CONTINUE 
100 CONTINUE 
RETURN 
END 

XDECK FFTY 

SUBROUTINE FFTY( SIGN, C0EF3 ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C FAST FOURIER TRANSFORM IN Y-DIRECTION X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

XCALL DATA9 
XCALL DATA* 

XCALL DAT21 

ISN=-SIGN 
DO 100 1=1, IMAX 
DO 110 J=1 , JMAX 
XR( J ) =FR(I , J ) 

XI ( J ) =FI (I , J ) 

110 CONTINUE 

CALL FFTCXR, XI, JMAX, ISN) 

IF (SIGN .LT. 0. ) GO TO 200 
DO 120 J = 1 , JMAX 
FR( I , J )=XR( J ) 

FICI, J)=XI( J) 

120 CONTINUE 
GO TO 100 

200 DO 130 J = 1 , JMAX 

FR( I , J ) =XR( J )XC0EF3 
FI ( I , J ) =XI ( J ) XC0EF3 
130 CONTINUE 
100 CONTINUE 
RETURN 
END 

XDECK FIX 

SUBROUTINE FIX( IM1 , I , IP1 , IMAX) 

IM1=I-1 

IP1=I+1 

IF(I .EQ. 1) IM1=IMAX 
IF. (I .EQ. IMAX) IP1 = 1 
RETURN 
END 

XDECK INICON 
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SUBROUTINE INIC0N( C , COF, DT , UR, VR, UR, UI , VI , WI , LI , 12 , 13 ) 

REAL NDIV , N12 > NSQR 

DIMENSION UR(L1,L2,L3),VR(L1,L2,L3) , WR( LI , L2 , L 3 ) 

DIMENSION UI(L1,L2,L3),VI(L1,L2,L3),WI(L1,L2,L3> 
c *** * X X X X X X X X X X XX X X X X X * X XX XX XXX X X X X X X XXXXXXXXX X X X X XX# X XXX X XXX XXX X XXXXXX X 
C THIS SUBROUTINE INITIATES THE PROGRAM, FOR STARTING PROBLEM, THE INI-* 

C TIAL FIELD IS GENERATED. FOR CONTINUATION PROBLEM, THE DATA STORED * 

C ON TAPE AT THE END OF THE PREVIOUS RUN ARE READ IN. * 

CK*K*HKKKK* XXXXXXXXX XX XXXXXX XX X XX XXXX X XX XXXXX XXXX XX XXX XXXXXX XXX XXX XXXXXX 

C NSTART=1 STARTING FROM TIME STEP=0 

C -NSTART = 2 CONTINUED FROM PREVIOUS RUN 

C IMAX=MAXIMUM MESH NUMBER IN X-DIRECTION 

C — JMAX=MAXIMUM MESH NUMBER IN Y-DIRECTION 

C LMAX=MAXIMUM MESH NUMBER IN Z-DIRECTION 

C TSTART=STARTING TIME STEP 

C TEND-ENDING TIME STEP 

C DELTA= MESH SIZE 

C- DT = TIME STEP 

C C=M0DEL CONSTANT 

C -NAVG=DELTA(AVERAGING)/DELTA(MESH) 

C ANISO= R IN EQUATION (5.16B) 

C GAMMA = STRAIN RATE 

REAL NAVG 

COMMON/ LARGEl/ENC 102 A) >EN1(102A),WN(20A8) 

LEVEL 2 , UR , VR , MR , U I , V I # WI 
LEVEL 2 , EN, EN1 , WN 
*CALL DATA7 

READ A, DELTA, DT,C, NAVG, ANISO, UTM, GAMMA 

IMAX=16 

JMAX=16 

LMAX=16 

IJK=IMAX*JMAX*LMAX 
I J=IMAX* JMAX 
IMMl-IMAX-1 
IMP 1=IMAX+1 
NMESH-IMAX 

TDIV=1 . 0/ C IMAX* JMAX*LMAX) 

NHALF=NMESH/2 

NHP1=NHALF+1 

HALF=FLOAT<NHALF) 

NM1=NMESH-1 
RIS0=3 ./(3.+ANIS0) 

TEMP=ANISO/3. 

RANISO=SQRT (TEMP) 

CC = 1. 

TFAC=IMAX*JMAX*LMAX 
FAC=SQRT ( TFAC) 

COEF10=3. 1A15926535898/NHALF 

C0EF1 1=C0EF1 0 

C0EF12=3.I A 15926535898*2. 

C0NST=COEF10/ DELTA 
C0EF15=C0EF12*FAC 
PI1-COEF10 
PI2=PI1*2. 

C0F=1 . 0 
NCONT = 1 
DO 2 M=l,25 
Y9=RGEN(X9) 

2 CONTINUE 

C ENERGY SPECTRUM DATA 

C- 0.1 INTERVAL UP TO 1.0 THEN .5 INTERVALL UP TO 6.0 

C EN IS THE ENERGY SPECTRUM FOR THE ISOTROPIC PART. EN1 IS FOR THE 

C ANISOTROPIC PART. 

PRINT 1960 

WN( 1 ) = 0 . 1 ,-:v 

DO 1900 M=2 , 1 0 
1900 WN(M)=0.1+WN(M-1) 

DO 1950 M-ll, 2A 
1950 WN(M)=0.5+WN(M-1) 

1960 FORMAT ( /5X, * WAVE NUMBER*,/) 










PRINT A, <WN(M) ,M=1,2A) 

PRINT 2000 

2000 FORMAT ( /5X » *UNF I LTERED SPECTRUMK,/) 

DO 3 M-l ,|2A ,8 
M7=M+7 

READ A, ( EN ( MM ) , MM=M > M7 ) 

PRINT 702 , ( EN(MM) » MM=M,M7 ) 

3 CONTINUE, 

DO 503 M»1,2A,8 
M7=M+7 i. 

READ A, ( EN1 (MM) ,MM=M ,M7 ) 

PRINT 702 , (ENl(MM) ,MM=M,M7 ) 

503 CONTINUE 

A FORMAT (8E10.A) 

DELAVG=(DELTAXNAVG)X*2/12.0 
PAI=5.1AI5926535898 
DO 2100 M= 1 , 2A 
EXPF=EXP(-DELAVGXWN(M)*WN(M)) 

EN(M)=EN(M)*EXPF 
2100 EN1(M)=EN1(M)*EXPF 
PRINT 2200 

2200 FORMAT (/5X, X FILTER ED SPECTRUM*,/) 

PRINT 702, (EN(M) ,M=1,2A) 

PRINT 702, (ENl(M) ,M=1»2A> 

DO 5 L-l/LMAX 
DO 5 J=1,JMAX 
DO 5 1=1, IMAX 
URL I , J , L ) = 0 . 

VR( I , J , L ) = 0 . 

WR(I,J,L)=0. 

UI ( I , J , L ) = 0 . 

VICI, J,L)=0. 

Uil ( I , J , L ) ~0 . 

5 CONTINUE 

DO AO L=1,NHALF 
LL=L 
N3=L-1 
N3S=N3**2 
DO 30 J=1 , NM1 
JINDEX=J/NHALF 
J J=J+NHP1-J INDEX* JMAX 
N2=J-NHALF 
N2S=N2**2 
DO 20 1=1, NM1 
I INDEX = I/NHAL F 
II=I+NHP1-IINDEX*IMAX 
N1=I-NHALF 
N15=N1**2 
NSQR=N1S+N2S+N35 
IF (NSQR . LT. 0.1) GO TO 20 
WAVN=SQRT (NSQR) 

NDIV=1 . /WAVN 
N12=N1S+N2S 

IF ( N12 .LT. 0.1) NC0NT=2 

IF (IABS(N1) .EQ. NHALF . AND. IAB5CN2) .EQ. NHALF) NC0NT=2 

C- GET FOURIER AMPLITUDE OF THE INITIAL FIELD AD DESCRIBED IN SEC A. A 

X=CONST*WAVN 
NREG=X+1 . 

GO TO (310,315,315,315,315,315,315) NREG 
310 M=X/0.1 

YM=X-0 .1XM 
M1=M+1 

ED=EN(MI)-EN(M) 

ENERGY=EN(M)+ED*YM*10. 

EDA=EN1(M1)-EN1(M) 

EANIS0=EN1(M)+ EDAXYMX10 . 

GO TO 320 
315 M=(X-1. )*2. 

YM=X~1 . “0 . 5*M 
M=M+10 
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M1=M+1 

ED=EN(M1 )~EN(M) 

ENERGY* EN(M)+ED*YM*2 . 

EDA=EN1 (Ml )~EN1 (M) 

EANIS0=EN1(M)+EDA*YM*2. 

320 Q5 = ENERGYXRIS0/ ( C0EF15XXXX2 ) 

QN=SQRT(QS) 

QSA=EANIS0XRIS0/(C0EF15XXX*2) 

QNA*$QRT( QSA ) 

CHANGE WAVE NUMBER VECTOR TO SATISFY NUMERICSl DIV FREE 

R1.R2 AND R3 ARE THE MODIFIED WAVE NUMBER 

IF( NCONT . EQ . 2 ) GO TO 340 

ARG1=PI1*N1 

ARG2=PI2XN1 

R1=ARG1/DEITA 

ARG1=PI1XN2 

ARG2=PI2*N2 

R2=ARG1/DELTA 

ARG1=PI1*N3 

ARG2=PI2*N3 

R3=ARG1/DELTA 

R1S=R1**2 

R2S=R2**2 

R3$=R3XX2 

R12S=R1S+R2S 

RSQ=R12S+R3S 

IF (NCONT. EQ. 2) GO TO 340 
R12=SQRT (R12S) 

R12DIV=I ./R12 

R12=SQRT(R12S) 

R12DIV=1 ./R12 
RDIV=1 ./SQRT(RSQ) 

GET A & B VECTOR 

FIRST CHOOSE RANDOM PHI 
340 CONTINUE 
YY=RGEN (XX) 

PHI=YY*C0EF12 

CPHI-COS(PHI) 

SPHI-SIN(PHI) 

I F (NCONT . EQ . 2 )GO TO 11 

A1=(-R2*CPHI+R1XR3*RDIVXSPHI)*R12DIV 

A2=(R1XCPHI+R2*R3*RDIV*SPHI)XR12DIV 

A3=-R12*RDIV*SPHI 

CALI RANDOM PHI 

Y2=RGEN(X2) 

PHI=Y2*C0EF12 

CPHI=COS(PHI) 

SPHI=S IN ( PHI ) 

B1=(-R2*CPHI+R1XR3*RDIVXSPHI)XR12DIV 
B2=(R1XCPHI+R2XR3*RDIV*SPHI)*R12DIV 
B3=-R12*RDIV*SPHI 
GO TO 12 

11 CONTINUE 
INDEX=(YY+0.25) *4 

PHI=0 .7853982X(2*INDEX-1) 

A 1 = S I N ( P H I ) 

A2=C0S ( PHI ) 

A 3 = 0 . 

Y1=RGEN(X1 ) 

INDEX=(Y1+0.25)X4 

PHI=0.7853982X<2XINDEX-1) 

B1=SIN(PHI ) 

B2-C0S(PHI ) 

B3 = 0. 

NCONT=l 

12 CONTINUE 

DETERMINE A AND B IN EQUATION (4.6) 

RANDOM THETA 
Y3“RGEN(X3) 
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URAN=QNA*CAKRANISO 
WIAN=QNA#CB*RANISO 

WR( 1 1 * J J * LL )=WR( II * J J > LL )+WRANXWSIGN 
WI(II,JJ,LL)=WI(II,JJ,LL)+WIANXVSIGN 
20 CONTINUE 
30 CONTINUE 
AO CONTINUE 

NOW THE UPPER HALF OF THE K-SPACE HAS BEEN DETERMINED 
GET THE TRANSFORMED VELOCITY AT THE CONJUGATE POINTS 

CONJUGATE FORM 

N3=l TO 7, N1 * N2=-7 TO 7 
N3=L-1 -N3=LM 

N2=J-1 -N2=JM 

IMP2=IMAX+2 
DO A1 L=2,NHALF 
LM=L+IMP2-2XL 
DO A 1 J = 1 > JMAX 
M=( J+IMM1 )/IMPl 
JM=J+(IMP2-2XJ)XM 
DO A1 I =1 » IMAX 
M=( I+IMM1 )/IMPl 
IM=I+( IMP 2-2X1 )XM 
URIIM, JM*LM>= URCI,J,L) 

VRUM, JM*LM> = VRCI*J,L) 

WRCIM,JM,LM)= WRCI,J,L) 

UKIM* JM,LM)=-UI(I, J,L ) 

VI ( IM* JM* LM) =-VI( I * J * L ) 

WI ( IM* JM* LM)=-WI( I * J * L ) 

A1 CONTINUE 
C N 3 = 0 * N1 = 1 TO 7, N2=-7 TO 7 

DO A2 1=2 ,NHALF 
IM=I+C IMP2-2XI ) 

DO A2 J = 1 * JMAX 
M=(J+IMM1)/IMP1 
JM=J+(IMP2-2XJ)XM 
IF(J.EQ.NHPl) GO TO A2 
URCIM*JM,1 )= UR(I,J*1> 

VR(IM*JM,1 )= VRCI*J*1) 

WR(IM*JM*1 > = WR(I,J*1) 

UKIM*JM*1 )=-UI(I,J*l> 

VI ( IM* JM* 1 )=-VI(I,J,l> 

WI ( IM* JM* 1 )=-WI(I,J,l) 

A2 CONTINUE 
C N1=N3=0 

DO A3 J=2*NHALF 
JM= J + (IMP2-2XJ) 

UR( 1 * JM* 1 ) = UR ( 1 > J * 1 ) 

VR ( I , JM , 1 ) = VR C 1 > J * 1 ) 

WR( 1 * JM* 1 ) = WR (1* J * 1 ) 

UI(l,JM,l)=-UI(l, J,l) 

VI(1,JM,1)=-VI(1,J*1) 

WI( 1 » JM* 1 )=-Wl ( 1 * J » 1 ) 

A3 CONTINUE 

C X AND Y TRANSFORM 

SIGN=-1. 

DO 50 L = 1 * LMAX 

CALL MOVL EV( UR( 1 * 1 * L ) * FRC1 * I ) * I J ) 

CALL MOVLEV(UI(l,l,L),FI(l,l),IJ> 
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THETA=Y3XC0EF12 

CA=COSCTHETA) 

CB=SIN(THETA) 

UR( II * J J , LL J = QNMCA*A! 
VR(II* JJ*LL)=QNXCAXA2 
WR(II*JJ,LL)=QNXCAXA3 
UICII*JJ,LL)=QNXCBXB1 
VI(II*JJ*LL)=QNXCBXB2 
WKII, JJ,LL)=QNXCBXB3 
IF (N3 .NE. 0) GO TO 20 
WSIGN=ABS(A3)/A3 
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CALL FFTX< SIGN ) 

CALL FFTY ( S IGN , CC ) 

CALL MOVLEV ( FRC 1 , 1 ) , URC 1 , 1 , L ) » IJ ) 

CALL MOVLEVCFICl,l),UI(l,l,L),IJ) 

CALL MOVLEVC VR(1 , 1 > L) , FRC 1 > 1 ) , IJ ) 

CALL M0VLEV(VI(1,1,L),FIC1,1),IJ) 

CALL FFTXC SIGN) 

CALL FFTY ( SIGN > CC) 

CALL MOVLEV (FR(1,1),VR(1,1,L),IJ) 

CALL MOVLEV (FI(I,1),VI(1,1,L),IJ) 

CALL MOVLEV ( WR ( 1, 1, L),FR( 1,1), I J) 

CALL MOVLEVC WI C 1 , 1 » L ) , FI Cl, 1 > , I J ) 

CALL FFTXC SIGN ) 

CALL FFTYCSIGN,CC) 

CALL MOVLEV < FR C 1 , 1) , WRC 1 , 1 , L ) , I J > 

CALL MOVLEVC FI (1,1), WIC 1,1* L),IJ) 

50 CONTINUE 

C 2 TRANSFORM 

DO 51 1 = 1 , IMAX 
DO 52 L = 1 , LMAX 
DO 52 J=1 , JMAX 
FRC J> L ) =UR ( I , J » L ) 

FI ( J » L ) =UI C I » J , L) 

52 CONTINUE 

CALL FFTY C SIGN , 1 . 0 ) 

DO 53 L = 1 , LMAX 
DO 53 J = 1 , UMAX 
UR(I,J,L)=FRCJ,L) 

53 CONTINUE 

51 CONTINUE 

DO 54 1 = 1 > IMAX 
DO 55 L=1,LMAX 
DO 55 J=l, JMAX 
FRC J , L ) =VR ( I , J , L ) 

FICJ,L)=VI(I, J,L) 

55 CONTINUE 

CALL FFTY (SIGN >1.0) 

DO 56 L =1 > LMAX 
DO 56 J =1 , JMAX 
VR C I > J > L ) =FR( J , L ) 

56 CONTINUE 

54 CONTINUE 

DO 57 1=1, IMAX 
DO 53 L=1 , LMAX 
DO 58 J = l, JMAX 
FR ( J > L ) =WR ( I > J > L ) 

FI C J > L ) =WI ( I , J , L) 

58 CONTINUE 

CALL FFTY C SIGN , 1 . 0 ) 

DO 5 9 L = 1 , LMAX 
DO 59 J=1 , JMAX 
WRCI, J,L)=FR(J,L) 

59 CONTINUE 

57 CONTINUE 

THE INITIAL FIELD HAS BEEN GENERATED. THE FOLLOWING IS TO PRINT 

OUT INFORMATION ON THE GENERATED FIELD 
VELOCITIES ARE STORED IN UR. VR AND WR 
-----TURBULENT ENERGY CHECK 
TKU=0 . 

TKV=0 . 

TKW=0. 

DO 95 L=l, LMAX 
DO 95 J=1 , JMAX 
DO 95 1=1, IMAX 
TKU=TKU+UR(I, J,L)X*2 
TKV=TKV+VRCI, J,L)x *2 
TKW=TKW+WRCI, J,L)XK2 
95 CONTINUE 

tku=tku*tdiv 

TKV=TKV*TDIV 
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TKW=TKW#TDIV 

TKSUM=TKU+TKV+TKW 

TKUR=TKU/TKSUM 

TKVR-TKV/TKSUM 

TKWR=TKU/TKSUM 

PRINT 707 

PRINT 706 

PRINT 700,DT,DELTA,C,NAVG,ANISO,UTM, GAMMA 

PRINT 702, TKU, TKV, TKW, TK5UM 

PRINT 702, TKUR,TKVR , TKWR 

PRINT 706 

PRINT 601. 

UT0T=0. 

VTOT-O . 

WT0T=0. 

DO 120 1=1 , LMAX 
PRINT 710, L 
USUM = 0 . 

VSUM=0 . 

WSUM = 0 . 

DO 116 J = 1 , JMAX 
DO 116 1=1 , IMAX 
U5UM=USUM+UR(I, J,L) 

VSUM=VSUM+VR(I,J,L) 

WSUM=USUM+UR( I , J , L ) 

116 CONTINUE 

PRINT 7 02, C URC 1,10, L), 1 = 1, NHAIF) 

PRINT 702 , ( VR( 1 , 1 0 , L ) , I =1 , NHAL F) 

PRINT 702, (UR (I, 10, L), 1=1, NHALF) 

PRINT 702, USUM, VSUM,U5UM 
UTOT=UTOT+USUM 
VTOT=VTOT+VSUM 
UTGT=UTOT+WSUM 
120 CONTINUE 

PRINT 702.UTOT, VTOT,UTOT 

700 FORMAT ( IX ,XINITIAL CONDITION. DT=X1PE1 0 . 4 , X DELTA=X1PE10 . 4, 

1 x G=x , 0PF7 . 4 » 3X> X AVER AGING GRID = x,F4.1, X DELTAX , / , 18X, 

2XANIS0=X,E12.5,3X,XUTM=X, E12.5,3X,XGAMMA=X,E12.5) 

702 FORMAT ( 1P8E15 . 7 ) 

705 FORMAT ( IX , XCONTINUED AT TIME STEP=x , 14 ,/ ,/ ) 

706 FORMAT ( 1H0 ,130HXXXXXXXXKXXXXXXXXXXXK#XXKXXXXXXXXXKXXXXXXXXXX*XXXXX 
-XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
-XXXXXXXXXXXX ) 

707 FORMAT (1H1) 

710 F0RMAT(1X,XPLANE=X,I3) 

711 FORMAT ( IX, X INITIAL CONDITIONX ,/lX, XDT = X1PE1 0 . 4 , X DELTA = X1 PEI 0 . 4 
- , X C=X , F7 . 4 , X U0=X1PE10 .4,/) 

601 FORMAT ( IX, XUM, VM, UMX ) 

RETURN 

END 

XDECK INIFILT 

SUBROUTINE INI FI LT ( U , N1 , N2 , N3 ) 

DIMENSION U(N1,N2,N3> 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C THIS SUBROUTINE IS USED TO FILTER THE INITIAL FIELD UITH A UIDE 
C FILTER ONLY IN THE Z-DIRECTION. 

CX XX XXXXXXXXXXXX XXX XXX XXXXXXXXXXXX XXX XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XX 

XCALL FLT 
XCALL DAT21 
XCALL DATA7 
XCALL DATA9 

LEVEL 2 , U 
LMAXM1 =LMAX~1 
DO 5 L = 1 , LMAX 

XR(L)=EXP(-FLOAT(L-1)XX2/8.0) 

5 CONTINUE 
AREA-0. 5XXR(1> 

DO 6 L=2 , LMAXM1 
AREA=AREA+XR( L ) 

6 CONTINUE 
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AREA=AREA+0.5XXR(LMAX) 

DO 7 L = 1 » LMAX 
XR(L)=XR(L)/AREA 

7 CONTINUE 

CALL FDCT (1.0) 

DO 8 L = 1 , LMAX 
FILT3(L)=XR(L) 

8 CONTINUE 

DO 1 J = 1 , JMAX 
DO 1 1 = 1 , IMAX 
DO 2 L = 1 , LMAX 
XR(L)=U(I, J,L) 

2 CONTINUE 

CALL FDCT < 1 . 0 ) 

DO 3 L = 1 , LMAX 
XR(L)=XR(L)*FILT3(L> 

3 CONTINUE 

CALL FDCT (-1.0) 

DO A L = 1 . LMAX 
UCI, J,L)=XR(L> 

A CONTINUE 
I CONTINUE 
RETURN 
END 

XDECK INVERS 

SUBROUTINE INVERS ( G, PM, HM, IC , N1 , N2 , N3 ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C INVERS IS A POISSON SOLVER. IC = 3 IT EXPANDS THE VARIABLE G 

C IN COSINE SERIES IN THE Z-DIRECTION OTHERWISE IT EXPANDS G 

C IN SINE SERIES IN THE Z-DIRECTION. IN THE OTHER TWO DIRECTIONS 
C FOURIER SERIES ARE USED TO EXPAND G 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DIMENSION G(N1,N2,N3),PM(N1,N2,N3>,HM(N1,N2,N3> 

XCALL DATA9 

XCALL DAT21 . 

XCALL DATA7 
XCALL WV 

LEVEL 2 r G, PM 
IJ=N1XN2 

CC=1 ./( IMAXXJMAX) 

CxxxxxTRANSFER G TO HM 
DO 10 J=1,JMAX 
DO 10 I =1 f IMAX 
DO 20 L = 1 , LMAX 
XR(L)=-G(I, J,L) 

20 CONTINUE 

IF (IC .EQ. 3) GO TO 100 
CALL FDST(l.O) 

GO TO 200 

100 CALL FDCT (1.0) 

200 DO 30 L=1 , LMAX 
HM(I, J,L)=XR(L) 

30 CONTINUE 
10 CONTINUE 

DO A0 L = 1 , LMAX 

CALL MOVLEV (HM( 1 , 1 , L ) » FR( 1 , 1 ) » I J ) 

CALL FFTX(l.O) 

CALL FFTY(1. 0,1.0) 

DO 50 J=1 , JMAX 
DO 50 1=1, IMAX 

WAV=WAVEXS ( I ) +WAVEYS ( J )+WAVEZS ( L ) 

IF (ABS(WAV) .LT. 0.00001) GO TO 500 
■■■■■■■-■■ WAV = 1 ./WAV 

FR(I,J)=FR(I,J)XWAV 
FI ( I , J ) =FI ( I , J ) XWAV 
GO TO 50 
500 FR(I,J)=0. 

FI ( I , J ) =0 . 

50 CONTINUE 

CALL FFTY ( -1 . 0 , CC) 
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CALL FFTX(-l.O) 

CALL MOVLEV<FR(l,l),HM(l,l,L),IJ) 

40 CONTINUE 

DO 60 J = 1 » UMAX 
DO 60 1=1 » IMAX 
DO 70 L=1,LMAX 
XR( L ) =HM( I > J » L ) 

70 CONTINUE 

IF (IC .EQ. 3) GO TO 300 
CALL FDSTC-1.0) 

GO TO 400 

300 CALL FDCTC-1.0) 

400 DO 80 L=l,LMAX 
PM( I , J , L ) -XR( L ) 

80 CONTINUE 
60 CONTINUE 
RETURN 
END 

XDECK MEANINI 

SUBROUTINE MEANINI 
COMMON/NORM/ DELU, THETA 
xCALL DEL 
C 

COMMON/L ARGE4/01D( 16 ,16, 33 ),02D( 16, 16, 33), 030(16,16,33) 

LEVEL 2,01D,02D,03D 
XCALL BLANK 
XCALL DIM 
XCALL LARGE2 
XCALL LARGE3 
XCALL LARGE5 
XCALL DATA9 

c XXX X X X X X X X XX XXX XX XXX XX X XX X X X X X XX X XX XXX X X X X X XXX X X X X X X XXX X X * XX X X X XX XXXXX X 
C THIS ROUTINE CREATS THE MEAN INITIAL FIELD . INITIAL SPICKS * 

C ARE STORED IN GU THEN FILTERED TO THE CREAT THE GAUSIAN CORE . * 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DO 500 L=1,LMAX 
DO 500 J=1 , JMAX 
DO 500 1=1, IMAX 
GU<I>J,L)=0. 

500 CONTINUE 

DO 501 J=1,JMAX 
GUC 6 , J , 17 ) =20. 

GU( 1 1 , J » 17 ) =20 . 

501 CONTINUE 
PRINT 1110 

1110 FORMAT ( 1H1 f 5X , * INITIAL VORTEX AT PLANE 1*) 

PRINT 1115 »■<■(' C GU C 1 , 1 , L ) » 1 = 1 , IMAX) » L) , L = 1 » LMAX) 

1115 FORMAT (IX, 16 F8. 2, 13) 

CALL STFILT 

CALL SFILTER<GU,DUDX,N1,N2,N3) 

PRINT 1113 

1113 F0RMAT(1H1,5X,*FILTERED VORTEX AT PLANE 1*) 

PRINT 1 1 15 , ( ( (GUC I » l , L) » 1= 1,16 ) , L ) , L=1 > LMAX) 

DO 508 L = 1 . LMAX 
DUDXt 1,1,1) =0.0 
DO 508 J = 1 > JMAX 
DO 508 1=1, IMAX 

DUDXC 1 , 1 , L ) =DUDX( 1 , 1 , L ) +GUC I , J , L )/( IMAX* JMAX) 

508 CONTINUE 

DO 509 L = 1 » LMAX 
DO 509 J=1 . JMAX 
DO 509 1=1, IMAX 
GUCI, J,L')=DUDXC1,1,L) 

509 CONTINUE 

DO 502 1 = 1, LMAX- 
DO 502 J=1 , JMAX 
DO 502 I =1 > IMAX 
02( I > J , L ) =02( I , J > L ) + GU( I , J , L ) 

502 CONTINUE 

CALL INVER5(01fGU»DUDXf1»N1»N2fN3) 
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CALL INVERSC02,GV,DUDX,2,N1,N2,N3> 

CALL INVERSC03,GW,DUDX,3,N1,N2,N3> 

CALL CURL0CGU,GV,GW,U,V,W,N1,N2,N3) 

COMPUTE THE MEAN INITIAL VELOCITY FIELD AND NORMILIZE THE EQUATION 
WITH DELTA U(DELU) FOR THE VELOCITY SCALE AND THETA THE MOMENTUM 
THICKNESS FOR THE LENGTH SCALES. 

CDUM=1 ./(IMAXKJMAX) 

DO 506 L =1 , LMAX 
DUDX( 1 , 1 , L ) =0 . 

DO 506 J = 1 > JMAX 
DO 506 1=1, IMAX 

DUDX(1,1,L)=DUDXC1,1,L)+UCI, J,L)XCDUM 

506 CONTINUE 

DELU=DUDX(1 ,1,LMAX)-DUDXC 1,1,1) 

THETA=C 0 « 25-*C DUDXC 1 » 1 , 1 )/DELU )#*2) *0 . 5 
LMAXMI =LMAX-1 
DO 504 L=2, LMAXMI 

THETA=THETA + C 0 . 25-C DUDXC 1 >1,L)/DELU)#X2) 

504 CONTINUE 

THETA=THETA+ C 0. 25- C DUDXC 1,1, LMAX)/DELU)XX2) *0.5 

THETA=THETA*DELTAZ 

DELTAX=DELTAX/THETA 

DELTAY=DELTAY/tHETA 

DELTAZ=DELTAZ/THETA 

CALL STWV 

DO 505 L = 1 > LMAX 

DO 505 J = 1 » JMAX 

DO 505 I=1,IMAX 

UCI, J,L>=UCI,J,L)/DELU 

VCI, J,L)=VCI, J,L)/DELU 

W(I,J,L)-’«I(I>J»L)/DELU 

505 CONTINUE j 

CALL CU^LUCU, V,W, OI, 02,03, N1,N2,N3) 

PRINT 708 

708 FORMAT ( 1H0, 13 OHXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXKX 
I X##X X X X X X X X XX XXX X X KK X X X XXX X X X X X XXX X X X X X X XX X X X XX X XX XX XX X X X# X XXX X X X X 
2XXXXXXXXXXXXX ) 

DT=0.03125XDELU/THETA 
PRINT 1116, DELU, THETA, DT 

1116 FORMAT C IX, 1H#, * DELU=* , E15 . 7 , 10X, * THETA=x , E15 . 7 , 1 OX , X DT=X,E15.7, 
446X,1H#) 

PRINT 1I17,DELTAX,DELTAY,DELTAZ 

1117 FORMAT! IX, 1 HX , X DELTAX=X,E15.7,10X>* DELTAY=x , E15 . 7 , 1 OX , * DELTAZ=x 
*,E15.7,35X,1H*) 

PRINT 708 
CALL RNDINIC 
UDUM=0. 

DO 507 L = 1 > LMAX 
DO 507 J = 1 , JMAX 
DO 507 1=1, IMAX 

IFC ABSCUC I , J, L) ) . GT . UDUM) UDUM = ABSCUC I , J , L ) ) 

IF(ABS(V(I, J,L)) .GT.UDUM) UDUM=ABS( VC I , J , L ) ) 

IFCABSCWCI, J, L)) .GT.UDUM) UDUM=ABS ( W( I , J , L ) ) 

507 CONTINUE 
CUDUM-0 .30/ UDUM 
DO 510 L=1 , LMAX 
DO 510 j =1 , JMAX 
DO 510 1=1, IMAX 

U(I, J,L)=UCI, J,L)#CUDUM 
V(I,J,L)=VCI,J > LJXCUDUM 
WCI. J,L)=WCI,J,L)xCUDUM 
510 CONTINUE 

CALL CURLUCU,V,W,01D,02D,03D,N1,N2,N3) 

DO 512 L = 1 , LMAX 
DO 512 J=l, JMAX 
DO 512 1=1, IMAX 

01 ( I , J , L > =01 < I , J , L ) + 01 DC I , j , L > 

02 C l , J , L ) =02C I , J , L ) + 02DC I > J » L ) 

03C I , J , L ) =03C I , J , L )+03DC I , J , L ) 

512 CONTINUE 



CALL INVERS<01,GU,DUDX,1,N1,N2,N3) 

CALL INVERS<02,GV,DUDX,2,N1,N2,N3> 

CALL INVERS(03,GW,DUDX,3,N1,N2,N3) 

CALL CURL0(GU,GV,GW,U,V,W,N1,N2,N3) 

DO 513 L = 1 » LMAX 
DO 513 J=1,JMAX 
DO 513 1 = 1 , IMAX 
OIDCI, J,L>=0. 

02DCI, J , L ) =0 . 

03D( I » J » L ) = 0 . 

513 CONTINUE 
DUMM1 =0 . 

DUMM2=0 . ' 

DO 3333 1=1, IMAX 
DO 3333 J=1 , JMAX 
DO 3333 L = 1 , LMAX 

IF(01(I»J,L). GT. DUMM1 ) DUMM1 =0 1 C I , J , L ) 

IF(02(I,J,L).GT. DUMM2 ) DUMM2=02( I , J , L ) 

3333 CONTINUE 

PRINT 333A » DUMM1 > DUMM2 

3334 FORMAT ( IX, X DUMM1= X,E15.7,2X,X DUMM2= X,E15.7> 

RETURN 

END 

XDECK PARPLOT 

SUBROUTINE PARPLOT 

c XX XX XX X XX X X X X X X XXX XX X XX X XX X X X XX X X X X XX X XXX X X XXX X X X X X XXX X XXX X XX X XX X X XXX X X 
C THIS ROUTINE PLOTS THE PARTICLES TRACKS . 

C XMIN IS FIXED TO BE ZERO 

C ZMIN IS FIXED TO BE ZERO 

C XMAX IS FLOATING AND DEPENDS ON NUMBER OF MESHES USED AND DELTAX 

C ZMAX IS FLOATING AND DEPENDS ON NUMBER OF MESHES USED AND DELTAZ 

C SCX IS THE SCALING FACTOR TO ADJUST TO A PAGE LENGHT OF 8 INCHES 

C SCZ IS THE SCALING FACTOR TO ADJUST TO A PAGE LENGHT OF 8 INCHES 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XCALL DATA9 
XCALL DEL 
XCALL XL 

DATA LB/1HX/ 

DATA NL/1HZ/ 

XMIN=0. 

XMAX=( IMAX-1 )x DEL TAX 
ZMIN=0. 

ZMAX=( LMAX-l )X DELTAZ 
PPXMAX=13 . 0 
PPZMAX=1 0 . 

SCX=PPXMAX/XMAX 

SCZ=PPZMAX/ZMAX 

CALL LINAXSCO. , 0 . , PPXMAX , PPZMAX, . 1 , -1 , 1 0 , 1 , XMIN , XMAX, 3 , 4 , LB ) 

CALL LINAXSCO. , 0 PPXMAX, PPZMAX, . 1 , +1 , 1 0 , 1 , ZMIN , ZMAX , 3 , A , NL ) 

DO 1 N =1,140 
X=XPART(N)XSCX 
Z-ZP ART ( N ) XSCZ 
NC=NCHARCN) 

CALL SYMBOL(X,Z,0.1,0,0.,-NC) 

1 CONTINUE 

CALL PLOTCO. ,0. ,6) 

RETURN 

END 

XDECK PARTRAC 

SUBROUTINE PARTRACCNPART»DT) 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

C THIS SUBROUTINE COMPUTES THE PARTICLES TRACK OF A TWO DIM MEAN. X 

C IT USES LINEAR INTERPOLATION TO COMPUTR THE VELOCITIES BETWEEN X 

C THE MESHES. TIME ADVANCING IS A FIRST ORDER EULER METHOD. X 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XCALL LARGE2 

XCALL DATA9 

XCALL DEL 

XCALL LARGE3 

XCALL XL 
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RLX=<IMAX-1)#DELTAX 

RLY=(JMAX-1)XDELTAY 

RLZ=(LMAX-1)XDELTAZ 

DO 1 M=1 , NPART 

IX=XPART(MVDELTAX+1 

IY=YPART(M)/DELTAY+1 

LZ=ZPART(M)/DELTAZ+1 

IXP1=IX+1 

IYP1=IY+1 

LZPI=LZ+1 




POOR 


IF( IX. EQ . I MAX) IXP1=1 
I F ( I Y .EQ. UMAX ) IYP1=1 
IF(LZ.EQ.LMAX) LZP1=LZ 
CCX=CXPART(M)-(IX-1)XDELTAX)/DELTAX 
CCY=( YP ART(M) -( IY-1 ) *DELTAY)/DELTAY 
CCZ=(ZPART(M)-(LZ-1)XDELTAZ)/DELTAZ 
U1PART=UCIX,IY , LZ >+(U(IXPl,IY 
,LZ ) + ( V( IXPI , IY 
, LZ )+(W( IXP1 » IY 
* LZP1 ) + <U( IXP1 , IY 
,LZP1)+CV(IXP1,IY 
, LZP1)+CW<IXP1,IY 


V1PART=V<IX,IY 
W1PART=W(IX,IY 
U2PART=U(IX,IY 
V2PART = V ( IX, I Y 
W2PART=W( IX, IY 


U1PART-U1PART+(U2PART-U1PART) XCCZ 


LZ )-U( IX, IY 
LZ )-V<IX,IY 
LZ ) -W( IX , IY 
LZP1)-UUX, IY 
LZPD-VCIX, IY 
LZPI )-UK IX, IY 


, LZ >)*CCX 
, LZ ))XCCX 

, lz nxccx 

, LZPI) )*CCX 
,LZP1))*CCX 
, LZP1 ) )*CCX 


W2PART=W(IX,IY ,LZP1)+(W(IXP1,IY , LZPI )-W( IX, IY ,LZP1))*CCX 
UlPART=UlPART+( U2PART-U1PART )*CCZ 
V1PART=VIPART+(V2PART-V1PART)XCCZ 
W1 P ART =W1 PART +( W2PART-W1PART } *CCZ 

U2PART”U(IX,IYP1,LZ ) + <U( IXP1 , IYPI , LZ >-U(IX, IYP1 , LZ ))*CCX 

V2PART = V<IX,IYP1»LZ )+ ( V( IXP1 , IYPI , LZ )-V ( IX, IYPI , LZ ))*CCX 

W2PART=W(IX,IYP1,LZ )+(W( IXPI , IYPI , LZ )-WC IX, IYPI , LZ ))*CCX 

U3PART=U( IX, IYPI, LZPI )+(U(IXPl, IYPI, LZPI )-U( IX, IYPI, LZPI ))*CCX 

V3PART=V(IX, IYPI, LZPI )+CV< IXPI, IYPI, LZPI )-V<IX, IYPI, LZPI >)*CCX 

W3PART=W( IX, IYPI, LZPI )+<W< IXPI, IYPI, LZPI )-W( IX, IYPI, LZPI)) XCCX 

U2PART=U2PART+(U3PART-U2PART)*CCZ 

V2PART = V2PART+ ( V3PART-V2P ART ) *CCZ 

W2PART=W2PART+(W3PART-W2PART)XCCZ 

U1PART=U1PART+(U2PART-U1PART)XCCY 

V1PART=V1PART+(V2PART-V1PART)XCCY 

W1PART=W1PART+CW2PART-W1PART)XCCY 

XPARTCM)=XPART(M)+DT*U1PART 

YPART(M)=YPART<M)+DT*V1PART 

ZPART(M)=ZPART(M)+DTXW1PART 

IFCXPART(M).GT.RLX) GO TO 10 

GO TO 20 

10 XPART(M)=XPART(M)-RLX 

20 IFCXPART(M) ,LT. 0 . ) GO TO 30 
GO TO AO 

30 XPART(M) =XPART(M)+RLX 

AO IFCYPART(M) .GT.RLY) GO TO 70 
GO TO 80 

70 YPART (M) -YPART(f1)-RLY 

80 IFCYPART (M) . LT . 0 . ) GO TO 90 
GO TO 100 

90 YPART(M)-YPART(M)+RLY 
100 IF(ZPARTCM) .GT.RLZ) GO TO 50 
GO TO 60 

50 ZPART (M) =RLZ 

60 IFCZPART(M) . LT . 0 . ) ZPART(M)=0. 

1 CONTINUE 
RETURN 
END 

XDECK PARTIAL 

SUBROUTINE PARTIAL CM, U , N1 , N2 , N3 ) 

DIMENSION U(N1,N2,N3) 

XCALL DATA9 
XCALL BLANK ^ 

XCALL WV 
XCALL DATA? 

. 'LEVEL- 2 , U 
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XX XX XX 


IJ=N1XN2 

IF CM .EQ. 2) GO TO 30 

CXXXXXDERIVATIVE IN THE X-DIRECTIONXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
DO 10 L = 1 » LMAX 

CALL MOVLEV(U(l»l,L),FR(l»l),IJ) 

CALL FFTX(l.O) f „ 

DO 15 J = 1 , JMAX J 

DO 15 1 = 1, IMAX 
DUM=FI ( I , J ) 

FI ( I , J ) =WAVEX( I)XFR<I»J) 

FR( I , J ) =-WAVEX( I ) XDUM 
15 CONTINUE 

CALL FFTXt-1.0) 

CALL MOVL EV ( FR( 1 , 1 ) > DUDX( 1 > 1 , L ) » I J ) 

10 CONTINUE 
GO TO 300 

CXXXXXDERIVATIVE IN THE Y- DIRECTION xxxxxxxxxxxxxkxxxxxxxxxxxxxxxxx 
30 CONTINUE 

DO 35 L = 1 , LMAX 

CALL MOV LEVY UC 1,1,L),FR(1»1),IJ) 

DO 32 J = 1 , JMAX 
DO 32 1=1, IMAX 
FI ( I , J ) = 0 . 0 
32 CONTINUE 

CALL FFTY( 1 . 0 , 1 . 0 ) 

DO 40 J=1 , JMAX 
DO 40 1 = 1, IMAX 
DUM=FI ( I , J ) 

FI(I,J)=WAVEY( J)XFR(I, J) 

FR( I , J ) =-WAVEY ( J ) XDUM 
AO CONTINUE 

CALL,, FFTYC-l. 0,1.0) 

CALL MOVL EV( FR( 1 > 1 ) » DUDX( 1 , 1 , L ) » I J ) 

35 CONTINUE 
300 CONTINUE 
RETURN 
END 

XDECK RGEN 

IDENT RGEN - PSEUDO RANDOM NUMBER GENERATOR 
XXX FUNCTION RGEN CD) 

CALLED AS A FUNCTION WITH 1 ARGUMENT (WHICH IS IGNORED) 
RETURNS IN X6 A RANDOM NUMBER GENERATED BY MULTIPLYING 
1 OF 5 INTEGER CONSTANTS BY THE CORRESPONDING GENERATOR 
SEE BKY USERS HANDBOOK FOR REFERENCES 

SST 

X RGENCOM - USED TO STORE THE GENERATORS AND POINTER 



USE 

/RGENCOM/ 


GEN 

DATA 

1048015011D 

THE 5 GENERATORS 


DATA 

2236846573D 



DATA 

42167 930 93D 



DATA 

7792106907D 



DATA 

963019 197 7 D 


PTR f 

DATA 1 


POINTER TO CURRENT GENERATOR 


USE 

X 



ENTRY 

RGEN 



IF 

-DEF, FTN, 1 



ENTRY 

RGEN$ 



COMMENT 

‘ RANDOM NUMBER GENERATOR (MMODLEVEL#) 

NAME 

VFD 

42/4LRGEN , 18/ RGEN 


IF 

DEF, FTN 



ELSE 

2 


RGEN 

PS 


ENTRY / EXIT 

RGEN$ 

EQU 

RGEN 

SINCE ARG IS IGNORED 


SA1 

PTR 

GET POINTER 


SA3 

RGEN 

GET ENTRY POINT 


SB 1 

1 



SX7 

4 



SA4 

Xl+GEN-1 

GET GENERATOR 


SA5 

Xl+CON-1 

GET CONSTANT 
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RGEN1 


MX2 

60-35 


AX3 

30 


IXO 

X7-X1 

A - PTR 

SB7 

X3 

B7 = RETURN ADDRESS 

DX6 

XAXX5 

GEN X CON 

PL 

XO , RGEN1 

UNLESS PTR =5 

MX1 

0 

IF PTR WAS 5 

BSS 

0 


SA3 

EXP 


SX7 

Bl+Xl 

INCREMENT PTR 

SA7 

PTR 

STORE NEW POINTER 

BX7 

-X2XX6 

MASK LOW 35 BITS 

BX5 

X7+X3 

PUT IN EXPONENT 

SA7 

AA 

STORE NEW GENERATOR 

NX6 

X5 


<JP 

B7 

JUMP DIRECTLY BACK 

DATA 

1 3 1 0 7 5 D 

CONSTANTS TO MULTIPLY BY 

DATA 

16 38 A 3D 


DATA 

196611D 


DATA 

2 2 9-37 9 D 


DATA 

2 6 21 A 7 D 


DATA 

173ABSA8 



CON 


EXP 

END 

XDECK RNDINIC 

SUBROUTINE RNDINIC 

CXXXXXXXXXX XX XX XXXXXXXXXXXXXXX XXXXXXXXXXXXXXXXXXXXX XXX X XX XXXXXXXXXXXXXXX 
C THIS SUBROUTINE CREATS THE INITIAL RANDOM FIELD BY CALL INICON X 

C INICON WAS INITIALY WRITTEN BY KWAK,D. AND IS USED HERE AS A MEANX 

C TO CREAT A RANDOM INITIAL FIELD. THE ROUTINE WAS WRITTEN FOR EQUALX 

C MESH, AND HENCE THE COMPLICATIONS IN THIS ROUTINE TO TRANSFER THE X 

C FIELD TO THE MIDDLE OF THE BOX . X 

C*XXKXX#XXXXXXXXXXXXXXX*XXXXXXXXXXXXXXXXXXXXX*HXXXX*XXXXXXXXXXXXXXX#XKXX 
XCALL BLANK 
XCALL LARGE2 
XCALL LARGE3 

COMMON/L ARGEA/Ol ( 16. 16, 33), 02 (16, 16, 33), 03(16,16, 33) 

LEVEL 2,01,02,03 
XCALL DATA9 
XCALL DIM 

C0MM0N/DUM1/ UMC 16 , 16 , 16 ) , VM( 16 , 16 , 1 6 ) , WMC 16 , 16 , 16 ) 

C0MM0N/DUM2/ GUIU6 , 16 , 16 ) , GVI ( 16 , 16 , 16 ) , GWI ( 16 , 16 , 16 ) 

LEVEL 2 , UM, VM, WM, GUI , GVI , GWI 

CALL INICON(C,COF,DT,UM,VM,WM, GUI , GVI, GWI, 16, 16, 16) 

IMAX=N1 
JMAX-N2 
LMAX=N3 
DO 1 1*1,16 
DO 1 J = 1 , JMAX 
DO 1 1=1, IMAX 
Ll=25-L 
L2=17-L 

U( I , J > LI ) =UM( I , J , L2 ) 

V ( I , J , L 1 ) - VM ( I , J , L 2 ) 

W(I,J,L1) =WM ( I , J , L 2 ) 

1 CONTINUE 

DO 2 L = 1 , 1 3 
DO 2 J=1 , JMAX 
DO 2 1=1, IMAX 
U(I, J,L)=0 . 

V(I,J,L)=0. 

W(I,J,L)=0. 

2 CONTINUE 

DO 3 L=21 , LMAX 
DO 3 J = 1 , JMAX 
DO 3 1=1, IMAX 
U(I,J,L) = 0. 

V < I , J , L ) = 0 . 

W( I , J , L ) = 0 . 

3 CONTINUE 

CALL INIFI LT ( U , N1 , N2 , N3 ) 
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CALL INIFILT(V,N1,N2,N3) 

CALL INIFILT<W,N1,N2,N3) 

CALL CURL0(U,V,W,01,02,03,N1,N2,N3) 

DO 16 L=1,LMAX 
DO 16 J=1,JMAX 
DO 16 1=1 i IMAX 
U( I , J , L ) =01 ( I , J , L ) 

V(I, J,L)=02(I, J,L) 

WCI, J,L)=03(I, J,L) 

16 CONTINUE 

CALL CURLU(U,V,W,01,02,03,N1,N2,N3) 

CALL INVERSC01,GU,DUDX,1,N1,N2,N3> 

CALL INVERS(02,GV,DUDX,2,N1,N2,N3) 

CALL INVERS ( 03 , GW, DUDX , 3 , N1 , N2 > N3 ) 

CALL CURL0(GU,GV,GW,U,V,W,N1,N2,N3) 

RETURN 

END 

XDECK SINPART 

SUBROUTINE S INPARTCU , N1 , N2 , N3 > 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C THIS ROUTINE COMPUTES THE PARTIAL DERIVATIVE OF U IN THE Z- 

C DIRECTION BY EXPANDING IN FOURIER SINE SERIES. 

C THE PARTIAL IS STORED IN DUDX. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DIMENSION U(N1,N2,N3) 

XCALL BLANK 
XCALL WV 
XCALL DAT21 
XCALL DAT A 9 

LEVEL 2 , U 
DO 10 J = 1 » JMAX 
DO 10 1=1, IMAX 
DO 20 L = 1 , LMAX 
XRCL ) =U( I > J , L ) 

20 CONTINUE 
SIGN=U4> 

CALL FDST ( SIGN ) 

DO 30 L=1 , LMAX 
XR(L)=XR(L)XWAVEZ(L) 

30 CONTINUE 
SIGN=“1 . 0 
CALL FDCT (SIGN) 

DO «0 L = 1 , LMAX 
DUDX ( I , J> L ) =XRCL ) 

AO CONTINUE 
10 CONTINUE 
RETURN 
END 

XDECK SFILTER 

SUBROUTINE SFILTER (HR, HI, N1,N2,N3) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXX 

C SFILTER FILTERS HR BY EXPANDING IT IN A FOURIER SINE SERIES IN X 

C THE Z/DIRECTIQN AND FOURIER SERIES IN THE OTHER TWO DIRECTIONS. X 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DIMENSION HR(N1,N2,N3),HI(N1,N2,N3) 

XCALL FLT 
XCALL DATA9 
XCALL DATA7 
XCALL DAT21 

LEVEL 2, HR 

CC=1 . 0/( IMAXX JMAX) 

IJ=N1XN2 
DO 10 J=1 , JMAX 
DO 10 1=1, IMAX - 
DO 20 L=1 , LMAX 
XR( L ) =HR( I , J , L ) 

20 CONTINUE 

CALL FDSTC1.0) 

DO 30 L = 1 , LMAX 
HI ( I , J , L ) =XR ( L ) 
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X X X 


30 CONTINUE 
10 CONTINUE 

DO 40 1=1 , LMAX 

CALL MQVLEV(HI(1>1»L)>FRC1»1),IJ) 

CALL FFTX(l.O) 

CALL FFTYCl. 0,1.0) 

DO 50 J=1,JMAX 
DO 50 1=1 , IMAX 

FR(I,J)=FR(I,J)XFILT1(I)XFILT2(J)XFILT3(L) 

FI ( I » J ) =FI ( I, J JXFILTl ( I )XFILT2< J )XFILT3( L) 

50 CONTINUE 

CALL FFTX(-l.O) 

CALL FFTYt-l . 0,CC) 

CALL MOVLEV(FR(l,l),HI(l,l,L),IJ) 

40 CONTINUE 

DO 60 J = 1 , JMAX 
DO 60 1=1, IMAX 
DO 70 L =1 , LMAX 
XR(L)=HI(I, J,L) 

70 CONTINUE 

CALL FDST(-l.O) 

DO 80 L = 1 , LMAX TTY OF THE 

HR ( I , J , L ) =XR ( L ) 

80 CONTINUE PXGE IS 

60 CONTINUE -PTGINAL 

RETURN 
END 

XDECK SG5 

SUBROUTINE SG5 ( U , V , E , N1 , N2 , N3 ) 

DIMENSION U(N1,N2,N3),V(N1,N2,N3),E(N1,N2,N3) 

XCALL DEL 
XCALL LARGE3 
XCALL LARGE5 
XCALL BLANK 
XCALL DATA9 

LEVEL 2 , U , V > E 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C THE SGS MODEL IS COMPUTED IN THIS ROUTINE BY SECOND ORDER DIFF X 

C AND STORED IN GU,GV,GW X 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
CSGSX=1./(2.XDELTAX) 

CSGSY=1 ./(2.XDELTAY) 

CSGSZ=1 . / ( 2 . XDELTAZ) 

IJK=N1XN2XN3 

CALL M0VLEVCDUDX(1,1,1),E(1,1,1) , IJK) 

DO 210 L = 1 , LMAX 
LM1 =L- 1 
LP1=L+1 

IF (L .EQ. 1) LM1 =LP1 
IF (L .EQ. LMAX) LP1=LM1 
DO 210 J=1 , JMAX 
CALL FIX( JM1 , J , JP1 , JMAX) 

DO 210 1=1, IMAX 

CALL FIX(IM1,I,IP1, IMAX) 

U( I, J,L)=CE(I, JP1,L)X01( I, JP1,L)-E(I, JM1,L)X01(I, JM1,L))XCSGSY 
I -(E(IP1,J,L)X02(IP1,J,L)-E(IM1, J,L)X02(IM1, J,L))XCSGSX 
V(I, J,L) = (E(I, J,LP1)X0KI, J,LP1)-E(I,J,LM1)X01(I, J , LM1 ) ) XCSGSZ 
1 - ( E( IP1 , J , L ) X03( IP1 , J , 1 ) -E( IM1 , J , L )X03 ( IM1, J , L ) ) XCSGSX 
210 CONTINUE 

CALL PARTI AL(2,U,N1,N2,N3) 

CALL M0VLEV(DUDXC1,1,1),UC 1,1,1), IJK) 

CALL C0SPART(V,N1,N2,N3) 

DO 220 L=1 , LMAX 
DO 220 J = 1 , JMAX 
DO 220 1=1, IMAX 

GU C I , J , L ) =GU( I , J , L ) +UC I , J , L ) +DUDXC I , J , L ) 

220 CONTINUE 

DO 230 L=1 , LMAX 
LM1=L-1 

LP1=L+1 ■ 



rswr-amni s s an-snr- , r pri &ZZ2233SSS :S fcfeg ■ ifegsa 't 


IF (L .EQ. 1) LM1=LP1 
IF (L . EQ. LMAX) IPI=LM1 
DO 230 J=1,JMAX 
CALL FIX( JM1 , J » JP1 , JMAX ) 

DO 230 1 = 1 , IMAX 

CALL FIXdMl,I,IPl,IMAX) 

U(I, J,L) = (E( Ipi , J > L)K02 (I P 1 > J,L)-EdMl, J,L)X02CIM1,J,L))»CSGSX 
1 ~<Ed, JPl + LJXOld, JP1,L)-Ed, JMl,L)*01d, JMl,l))*CSGSY 
, V( h;!; L >!, <E<I ' J ' LPi) * 02(I ' J ' LP1) - E<1 ' J ’ LM1, *02d.a,i.Mi)>xcsGsz 

230 1 CONflHUE P1,L)K ° 3<I ' JP1,L) ~ E<1 ' JMl ' L))<03(1 ’ Jf11 ’ L) >><CSGSY 
CALL PARTIAL ( 1 , U,N1 ,N2 , N3) 

CALL MOVLEV(DUDX< l,l,l),Ud,l,l),IJK) 

CALL COSPART (V,N1,N2,N3) 

DO 240 L = 1 , LMAX 
DO 240 J = 1 , JMAX 
DO 240 1=1, IMAX 

£)(<I;id>=GV<I’ J *L ) +Ua,J,L>+DUDX<I,J.l) 

240 CONTINUE 

DO 250 L = 1 , LMAX 

lMl=L-l 

lPl=L+l 

IF (L .EQ. 1) LM1=LP1 
IF < L .EQ. LMAX) LP1 = LM1 
DO 250 J = 1 , JMAX 
CALL FIX(JM1, J, JP1, JMAX) 

DO 250 1=1, IMAX 

CALL FIX( IM1 , I , IP1 , IMAX) 

U(I, J,L) = (E(I pi , J ,L)K03d pl , U^J-EdMl, J,L)X03(IM1,J,L))XCSGSX 
1 .,7i E< . 1 .' J'LPUxOidfJ.lPn-Ed, J,LMl)*01d, J , LM1 ) )XCSGSZ 
i V< h ^ L> , = . ( £H' JP1 ' L) * 03(I ' Jpi ,L)- E d, JM!,L)*03d, JM1,L))XCSGSY 

250 l C0NflNUE ,LP1)XO2a ' J ' LP1> ‘ E<I ' J ' LM1))<O2<I ' J ' LN1))KCSGSZ 
CALL PARTIALd,U,Nl»N2,N3) 

CALL MOVLEV(DUDXCl,l,l),Ud,l,l),IJK) 

CALL PARTIAL(2,V,N1,N2,N3) 

DO 260 L = 1 , LMAX 
DO 260 J=1 , JMAX 
DO 260 1=1, IMAX 

260 CONTINUE' =G *^^ ^ » j * l ^ ( I , j , L 14 DUDX (I , J , L ) 

RETURN 

END 

*DECK STFILT 

SUBROUTINE STFILT 

CXXXXXXXXXXXXKXKXKKKkkkkkkxkk###*#*##***#*#**##*#^#^^,^^^^^^^^ 

£ lTii,^y25 OU I, I ,t! E JNITIALIZE the TRANSFORM OF the FILTER IN EACH 
C ?H lu^RoShNE FUTEr" R " “ S, ° RED IK F t>-Tl . FILT2. FIIT3. FOR USE 
CKKK*fcKXXKKKXX*KXKKXttKKMMKKXXMttKXNXMttKXKX*ttKXKKKXKKKXMMMKXKMKXXtfttXXMXXKX 

"CALL AVG 


xCALL FLT 
XCALL DATA7 
XCALL DAT21 
XCALL DATA9 
XCALL PR 

NHP1X= IMAX/ 2+1 
NHP1Y= JMAX/ 2+1 
1 "r NHP2X=NHP1X+1 
r NHP2Y=NHP1 Y+l 
LMAXM1=LMAX-1 
IF(CCF .NE. 0 . ) GO TO 400 

CXXXXXFIX THE TRANSFORM OF THE FILTER IN THE X-DIRECTION 
DO 100 J=1,JMAX 
DO 100 1=1 ,NHP1X 

100 cont:nue EXP< " 6 * <FL0 * T<I ' 1>/ * WG1), " J> 

DO 110 J=1 , JMAX 
DO 110 I=NHP2X» IMAX 
II=IMAX-I+2 


190 


S 


FR( I » J ) =FR( II > J ) 

110 CONTINUE 
AREA=0 .0 
DO 120 1 = 1 > IMAX 
AREA=AREA+FR( I > 1 ) 

120 CONTINUE 

DO 130 J = 1 » JMAX 
DO 130 1=1, IMAX 
FRCI, J)=FR(I, J)/AREA 


OF THl 

ORIGINAL PAGE IS POOR 


FI ( I , J ) =0 . 0 
130 CONTINUE 

CALL FFTX(l.O) 

DO 140 1=1, IMAX 
FILT I C I ) =FR( 1,1) 

140 CONTINUE 

CXXXXXFIX THE TRANSFORM OF THE FILTER IN THE Y-DIRECTION 


DO 200 J = 1 , NHP1Y 
DO 200 1=1, IMAX 

FR( I > J ) =EXP( -6 . x( FLOAT ( J-l )/AVG 2 )xx 2 ) 


200 CONTINUE 

DO 210 J =NHP2Y , JMAX 
DO 210 1=1, IMAX 
J J = JMAX- J + 2 
FR(I,J)=FR(I,JJ) 

210 CONTINUE 
AREA=0 .0 
DO 220 J = 1 , JMAX 
AREA=AREA+FR(1, J) 

220 CONTINUE 

DO 230 J = 1 , JMAX 
DO 230 1=1, IMAX 
FRCI, J)=FR(I, Jl/AREA 
FI ( I , J ) =0 . 0 
230 CONTINUE 

CALL FFTYCl. 0,1.0) 

DO 240 J=1 , JMAX 
FILT2C J ) =FR( 1 , J ) 

240 CONTINUE 

CXXXXXFIX THE TRANSFORM OF THE FILTER IN THE Z-DIRECTION 
DO 300 L=1 , LMAX 

XRCL)=EXP(-6.X(FL0AT(L-1)/AVG3)XX2) 

300 CONTINUE 

AREA = 0 . 5XXRC 1 ) 

DO 310 L = 2 , LMAXM1 


AREA=AREA+XR( L ) 

310 CONTINUE 

AREA=AREA+0 .5XXRCLMAX) 

DO 320 L=1 , LMAX 
XR( L ) =XR( L )/AREA 
320 CONTINUE 

CALL FDCTC1.0) 

DO 330 L=1 , LMAX 
FILT3 ( L ) =XR( L ) 

330 CONTINUE 
‘ FILT 1 ( NHP1X) =0 . 

FILT2(NHP1Y)=0 . 

FILT3C LMAX) = 0 . 

GO TP 410 

400 IFCCCF .NE. 1.0) GO TO 410 
MC=( LMAX-1 )X2/3 
DO 7 L=1 , LMAX 
FILT3(L)=0. 

7 CONTINUE 

DO 8 L = 1 , MC 
iFI LT3 C L ) = 1 . 0 

8 CONTINUE 
MC= JMAX/3+1 
DO 9 J=1 , JMAX 


FILT2C J)=0. 
9 CONTINUE 


DO 10 J=1,MC 
FILT2CJ)=1.0 
JJ=JMAX-J+1 
FIIT2(JJ) = 1.0 ,( 

10 CONTINUE 
MC= JMAX-MC+1 
FILT2(MC)-0. 

MC=IMAX/3+l 

DO 11 1=1. IMAX 
FILT1 ( I ) = 0 . 

11 CONTINUE 

DO 12 1=1, MC 
FILTlCI)=i.O 
I I=IMAX+1-I 
FILT1 ( II ) = 1 . 0 

12 CONTINUE 
MC=IMAX-MC+1 
FILT1(MC)=0 . 

410 PRINT 1119, CCF 
1119 FORMAT C 1H0 , X CCF = x , 1PE15 . 7 ) 
IF(CCPF .NE. 1.0) GO TO 340 
PRINT 1116, (FILTl(I), 1=1, IMAX) 
PRINT 1117,(FILT2(J),J=1, JMAX) 
PRINT 1118, (FILT3(L), 1=1, LMAX) 

1116 FORMAT ( IX , x FILT1 =*,1PE15.7) 

1117 FORMAT ( IX, X FILT2 =X,1PE15.7> 

1118 FORMAT ( IX, X FILT3 =X,1PE15.7) 
340 CONTINUE 

RETURN 

END 

XDECK STPART 

SUBROUTINE STPART( NPART ) 

XCALL LARGE5 
XCALL XL 
XCALL DAT A 9 
XCALL DEL 

YPART ( 1 ) = 0 . 

YPART ( 17 ) =0 . 

XPART(1)=(6-1)XDELTAX 

XPART(17)=(11-1)XDELTAX 

ZPART(1)=<17.-1. )XDELTAZ 

ZPART(17)=(17.-1.)XDELTAZ 

NCHAR ( 1 ) = 1 

NCHARC 17 ) =2 

YPART C 33) =0 . 

YPART( 4 9 ) = 0 . 

YPART ( 65 ) =0 . 

YPART (81 ) = 0 . 

YPARTC 97 ) = 0 . 

YPART ( 1 13 ) =0 . 

YPART < 129 ) = 0 . 

YPART ( 145 ) = 0 . 

XPART(33)=XPART(l)+0 .5XDELTAX 
XPART<49)=XPART(1)-0.5XDELTAX 
XPART(65) =XPART ( 1 ) 
XPART(81)=XPART(1> 

ZPARTC33) =ZPART ( 1 ) 

ZPART(49) =ZPART ( 1 ) 
ZPART(65)=ZPART(l)+0 .5XDELTAZ 
ZPART(81)=ZPART(l)-0 .5XDELTAZ 
XPART(97)=XPART(17)+0.5XDELTAX 
ZPART(97)=ZPART<17) 

XPART( 113 ) =XPART ( 1 7 )-0 . 5x DEL TAX 
ZPART < 113 ) =ZPART ( 17 ) 
XPART(129)=XPART(17) 
ZPART(129)=ZPART(17)+0.5XDELTAZ 
XPART C145)=XPARTC17) 
ZPART(145)=ZPART(17)“0.5XDELTAZ 
NCHARC 33>=3 
NCHARC 49) =4 


192 





poor 


NCHARC 65 ) =5 
NCHARC 81 ) =6 
NCHAR( 97)=7 
NCHARC 1 1 3 ) =8 
NCHAR(129)=9 
NCHARC 145 ) =10 
NPART-160 
DELX'-O . 

DELZ=0. 

N = 0 

DO 1 M=l,10 
N=N+1 

DO 2 J=2, JMAX 
N=N+1 

NCHARCN) =M 

IX=XPARTCN-1)/DELTAX+1 

LZ=ZPART(N-1)/DELTAZ+1 

IXP1=IX+1 

LZP1=LZ+1 

CCX = C XPARTC N-l >- ( IX-1 >*DELTAX)/DELTAX 
CCZ=CZPARTCN-1)-CLZ-1)*DELTAZ)/DELTAZ 
01P1=01CIX,J,LZ)+C01CIXP1, J,LZ)-01CIX, J,LZ))*CCX 
02P1=02CIX, J,LZ)+C02CIXP1, J,LZ)-02CIX, J,LZ))*CCX 
03P1=03CIX, J,LZ)+C03CIXP1, J,LZ)-03CIX, J, LZ))*CCX 
01P2=01CIX, J,LZP1)+C01CIXP1,J,LZP1)-01CIX, J,LZP1) >*CCX 
02P2=02CIX,J,LZP1)+C02CIXP1,J,LZP1>-02CIX,J,LZP1))*CCX 
03P2=03CIX, J,LZP1)+C03CIXP1,J,LZP1)-03CIX, J,LZP1))*CCX 
01P1=01P1+C 01P1-01P2 ) *CCZ 
02P1=02P1+C02P1-02P2)*CCZ 
03P1 = 03P1+C03P1-03P2 ) *CCZ 
DELX=01P1*DELTAY/02P1 
DELZ=03P1*DELTAY/02P1 
XPART ( N ) =XPART CN-ll+DELX 
YPARTCN)=YPARTCN-1)+DELTAY 
ZPARTCN)=ZPART(N-1)+DELZ 
2 CONTINUE 
1 CONTINUE 
RETURN 
END 

*DECK STREAD 

SUBROUTINE STREAD 

c * *XX* XXX XXX * X X X K X X* XXX X X* XX X X XXXXXX X XX X X XXX# * XX X X XXX X X X XX XX X X X X X X XX X X X X 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS SUBROUTINE READ THE INPUT PARAMETERS 
IMAX=NUMBER OF GRID POINTS IN THE X-DIRECTION 


JMAX=NUMBER OF 
LMAX=NUMBER OF 
AVG2=FILTERING 
AVG3=FILTERING 
AVG1=FILTERING 
DELTAX= MESH SIZE 
DELTAY= MESH SIZE 
DELTAZ= MESH SIZE 
Nl= ARRAY SIZE IN 
N2= ARRAY SIZE IN 
N3 = ARRAY SIZE IN 
CCFW= 1 IF PRINT 
CCPF= 1 IF PRINT 


GRID POINTS IN THE Y-DIRECTION 
GRID POINTS IN THE Z-DIRECTION 
WIDTH IN THE Y-DIRECTION 
WIDTH IN THE Z-DIRECTION 
WIDTH IN THE X-DIRECTION 
IN THE X-DIRECTION 
IN THE Y-DIRECTION 
IN THE Z-DIRECTION 
THE X-DIRECTION 
THE Y-DIRECTION 
THE Z-DIRECTION 

OUT OF WAVE IS WANTED , OTHERWIZE 
OUT OF FILT IS WANTED, OTHERWIZE 


NO 

NO 


PRINT 

PRINT 


OUT 

OUT 


CCPD= 1 IF PRINT OUT OF LINE AVERAGE OF U-COMPONENT 

LINE AVERAGE OF W-COMPONENT AND ENSEMBLE AVERAGE PERTURBATIONS* 
IS REQUIRED , OTHERWIZE NO PRINT OUT . * 

C**************** ****************************************** ************* 
INTEGER TSTART , TEND 
*CALL DATA 9 

COMMON/TIM/ TSTART, TEND 
XCALL DEL , 

*CALL DIM 
*CALL AVG 
*CALL PR 

READ 703 , IMAX, JMAX, t MAX, TSTART , TEND 
READ 704, DEL TAX, DELTAY, DELTAZ 


READ 704,AVG1,AVG2,AVG3,CCF 
READ 703,N1,N2,N3 
READ 704,CCPW,CCPF,CCPD 
PRINT 708 

PRINT 705, IMAX, JMAX, LMAX, TSTART, TEND 
PRINT 706,DELTAX,DELTAY,DELTAZ 
PRINT 7 0 7 , AVG1 , AVG2 . AVG3 
PRINT 7 0 9 , N 1 , N 2 , N 3 
PRINT 708 
703 FORMAT (1015) 

7 04 FORMAT ! 4E1 0 . 4 ) 

705 FORMAT! IX, X IMAX=X , 15 , 5X, X JMAX=X, 15, 5X, X LMAX=X , 15 , 5X , X TSART=XI5 
+ ,5X,x TEND=X,I5,52X,1HX) 

706 FORMAT ( IX, X DELTAX=X , 1PE10 . 4 , 5X, X DELTAY=X , 1 E10.4,5X,x DELTAZ=X,1 
+ PE10 . 4 , 64X, 1HX) 

707 FORMAT ! IX , X AVG1X=X, 1PEI0 . A, 5X, X AVG 2 =x,l E10.4,5X,X AVG3=X,1 
+PE10.4.64X, 1HX) 

708 FORMAT! 1 H 0 , 130 HxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxXxxxxxxxxxx 

lXXXXXXXXXKXXXKXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXX 
2 ^ K K ^OOOOOOOOOC ) 

709 FORMAT! IX, X N1=X,I5,5X,X N2=*,I5,5X,X N3=K, 15, 5X, 70X, 1H*) 

RETURN 

END 

XDECK STWV 

SUBROUTINE STWV 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


C STWV SETS THE WAVE NUMBERS FOR A GIVEN MESH SIZE DELTA AND X 
C NUMBER OF MESH POINTS NMAX .THIS ROUTINE MUST BE CALLED X 
C TO INITIALIZE THE WAVE NUMBERS FOR THE PARTIAL ROUTINES AND X 
C INVERS ROUTINE . X 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

XCALL WV 
XCALL DATA9 
XCALL DEL 

XCA I I pc 

PAI=3. 1415926535898 

CX=2. 0XPAI/ ! FLOAT !IMAX)XDELTAX) 

CY=2. 0XPAI/! FLOAT! JMAX1XDELTAY) 

CZ=PAI/! FLOAT! LMAX-DXDELTAZ) 

C2X=CX/ FLOAT ! I MAX) 

C2Y=CY/FLOAT! JMAX) 

NHPlX=IMAX/2+l 
NHP1Y= JMAX/ 2+1 
DO 100 L=1 , LMAX 
WAVEZ! L ) =CZX FLOAT ! L~1 ) 

WAVEZS!U=-WAVEZ! L)XX2 

100 CONTINUE 

DO 101 J = 1 , JMAX 

MM=J/NHP1Y 

M=MMXJMAX+1 

WA V EY !J)=C2YX FLOAT !J-M) 

WAVEYS! J)=-!CY*FL0AT!J-M))XX2 

101 CONTINUE 

DO 102 1=1 , IMAX 

MM=I/NHP1X 

M=MMXIMAX+1 

WAV EX! I ) =C2XX FLOAT ! I-M) 

WAVEXS! I ) =- ! CXX FLOAT! I-M) )XX2 

102 CONTINUE 
WAVEXC NHP IX) = 0 . 

WAVEY! NHP1Y) =0 . 

WAVEXS ! NHP IX) =0 . 

WAVEYS!NHP1Y)=0. 

WAVEZ! LMAX ) =0 . 

WAVEZS !LMAX)=0. 

IF! CCPW .NE. 1) GO TO 104 

PRINT 1000 , !WAVEX!I) .WAVEXS! I), 1=1, IMAX) 

PRINT 1001, ! WAVEY! J), WAVEYS! J),J=1, JMAX) 

PRINT 1002, !WAVEZ!L),WAVEZS!L),L=1, LMAX) 


104 CONTINUE 

1 000 FORMAT < IX, X WAVEX =* , 1PE15 . 7 , 5X , X WAVEXS =X,1PE15.7) 

1001 FORMATS IX# X WAVEY = X , 1PE15 . 7 , 5X , X WAVEYS =X,1PE15.7) 

1002 FORMAT ( IX# X WAVEZ =X , 1PE15 . 7 , 5X , X WAVEZS =X,1PE15.7) 
RETURN 

END 






